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SUMMARY 


A two-dimensional advanced panel far-field potential flow model of the 
undistorted, interacting wakes of multiple lifting surfaces has been developed 
which allows the determination of the spanwise bound circulation distribution 
required for minimum induced drag. This model has been implemented in a 
FORTRAN computer program, the use of which is documented in this report. 

The nonplanar wakes are broken up into variable sized, flat panels, as 
chosen by the user. The wake vortex sheet strength is assumed to vary linearly 
over each of these panels, resulting in a quadratic variation of bound circu- 
lation. Panels are infinite in the streamwise direction. The theory is briefly 
summarized herein; sample results are given for multiple, nonplanar, lifting 
surfaces, and the use of the computer program is detailed in the appendixes. 

INTRODUCTION 


With the current resurgence of interest in utilization of unconventional 
aircraft concepts for future transport aircraft to provide reductions in drag, 
increases in fuel efficiency, and lower operating costs, there is a need for 
accurate estimations of the induced drag for nonplanar configurations. Examples 
of these novel configurations include wings fitted with end plates or winglets 
(ref. 1), the tandem wing (ref. 2), and the joined wing (ref. 3). The current 
far-field theoretical model allows very accurate calculation of induced drag 
for multiple nonplanar aerodynamic surfaces, allowing investigation of the 


drag-reduction potential of such nonplanar aircraft concepts. Further, the 
bound circulation and wake vortex strength distributions necessary to achieve 
this minimum induced drag are computed. The bound circulation output may then 
be used to determine the aerodynamic surface camber shapes required to achieve 
this minimum drag. 

The theoretical wake model has been described in detail in reference 4. 

It assumes the vortex sheet strength to vary piecewise linearly on a number 
of flat wake panels. Wake rollup is neglected. Analytical expressions are 
developed for the normal velocities induced by each wake panel at any point 
on the wake using the Biot-Savart Law (ref. 5). The wake vortex strength is 
integrated spanwise to compute bound circulation, and the product of local 
bound circulation with the total induced normal velocity is analytically 
integrated to obtain the induced drag. To determine the wake vortex distri- 
bution required for minimum drag, two theoretical methods are used: Munk's 

criterion (ref. 6), and a direct optimization technique. The Munk criterion 
technique is computationally more efficient, since only the induced velocity 
expressions are utilized. This technique is similar to the theory developed 
in reference 7. The direct optimization technique is necessary for determi- 
nation of induced drag for relative optimum configurations which might have 
additional constraints on bending moment or pitching moment. Analytical 

expressions are developed for the derivatives of C , C , C , and C with 

D L m B 

respect to the unknown values of the wake vortex sheet strengths at the 
corners of each vortex sheet panel, as described in reference 8. This wake 
model using the direct optimization technique has been implemented in a vortex 
lattice wing design computer code (ref. 9), as described in reference 8. 
Comparisons between results of the original design code (ref. 9) and modified 
code with the current wake model (ref. 8) are given in references 8 and 10. 

This report details the computer program which was written to implement 
the theoretical wake model of reference 4. The theory is briefly summarized. 
Use of the program and sample input and output data are given in the appendixes 
the code is briefly described (Appendix A) , input data preparation is explained 
(Appendix B) , output data is described (Appendix C) , sample input and output 
data are given (Appendix D) , and a listing of the computer program is given 
(Appendix E) . 
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SYMBOLS 


A, A. . 
ij 


matrix of influence coefficients in induced drag 
[eq. (11)] 



c 



C 


B 


C 


D 



C 


L 


C 

m 

d 


G. 

1 




i 


h 



h! . 
ij 


integrals appearing in normal wash expression 

[eq. (3)] 

constant appearing in integrals in Appendix A 
wing span 

constant appearing in integrals in Appendix A 
pressure coefficient 
difference in pressure coefficient 
wing root bending moment coefficient 
induced drag coefficient 

induced drag coefficient on wake panel i due to induced 
velocity of panel j and its image [eq. (6)] 

lift coefficient 

pitching moment coefficient 

constant appearing in integrals in Appendix A 

variables containing unknown wake vortex sheet strengths, 
appearing in drag coefficient equation [eq. (6)] 

vertical separation of diamond wing roots (figs. 6, 7) 

distance between influenced point on panel i and influencing 
point on panel j (fig. 1) 

distance between influenced point on panel i and influencing 
point on image of panel j (fig. 1) 


I I. I 2. . ’ 

1.3 1,3 

I 3. I 4. 

1.3 i,3 



influence coefficient integrals appearing in drag coefficient 
equation [eq. (6)] 
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N 

R 

R. . 
J-J 

R'. . 

ij 



s 

T 

U 

w 

n 
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w 

o 

X 

Y 
Z 

Y 

r 

7 

r 

o 

p 

X 

°±3 
<P ' 


induced drag efficiency factor, defined as ratio of planar 
wing induced drag to that of nonplanar configuration 

vertical fence height (fig. 2) 

integer appearing in integrals in Appendix A 

number of wake panels on one-half of total configuration 

constant appearing in integrals in Appendix A 

projection of distance h.. onto the plane of influenced 
panel i (fig. 1) 

projection of distance h! . onto the plane of influenced 
panel i (fig. 1) 

reference wing area 

local wake panel coordinate 

wake panel semiwidth 

constant appearing in integrals in Appendix A 
free-stream velocity 
normal wash velocity 

normal wash velocity induced at point 4=4^ on panel i 
to panel j and its image [eq. (2)] 

constant appearing in Munk’s criterion normal wash velocity 
expressions [eqs. (7) and (8)] 

streamwise coordinate 

spanwise coordinate 

vertical coordinate, positive down 

wake trailing vortex sheet strength 

bound circulation 

average bound circulation 

bound circulation at outboard edge of wake panel [eq. (4)] 

nondimensional spanwise coordinate 

Lagrange multiplier in equations (10) and (12) 

angle between y-axis and orientation of h„ (fig. 1) 

•dihedral angle 


due 
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Subscripts : 
i 

j 

n 

o 


influenced point 
influencing point 
normal component 
value at 4 = -s 


Superscripts : 

1 image panel quantity- 

average value 

THEORETICAL DEVELOPMENT 


The wakes are broken into N finite, flat panels. The panels are 
numbered sequentially from tip to root of each planform. On each panel, 
the wake vortex sheet strength is assumed to vary linearly, as 




Vi + Y i 




( 1 ) 


The y value equals the vortex sheet strength at the junction between 
segments j and j-1 (see the inset on fig. 1). Figure 1 also shows the wake 
geometry notation used in the current theory for a wing-winglet configuration. 


The induced velocity normal to point 4^ on wake segment i, due to 
wake segment j and its image, may be written, using the Biot-Savart Law, as 



( 2 ) 


The R„ and h„ are distances as shown in figure 1. This expression is 
evaluated (ref. 4) as 
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w . = Y j+1 + ^ 
n,J o 


1 ( A 1.. +A 2..) +Iiti r Jl ( A 3.. +A 4..) 
IJ ir ' IJ IJ 


( 3 ) 


where the , etc. , are determined by the wake geometry. Next, the y 

ij 

distribution is integrated spanwise beginning at the tip of the current 
planform to obtain the bound circulation as 

' Y • , . + Y,\ 


r( 4 i ) ■ r of s i) + T (n +1 + a) + ( -" 2 - i )" 1 

. ( y i+l - T i \ *1 

\ 2 /2s. 


(4) 


where r ( -s . ) is the value of r at 4 . = -s . , which is a known linear 
o\ 1 / x i 

function of the y values. 

The drag induced on segment i by segment j and its image is, in 
coefficient form. 


+s. F(4.) /w . (4.X 

/ 1 -IT- (- 




D,ij S 


ref -s. 

l 


(5) 


This equation has been evaluated analytically as described in reference 4 
using the MACSYMA symbolic manipulation language (ref. 11). The result is 
an expression which is quadratic in the unknown y. values: 


C . . = — — - — ("g.G.I, + G.G.I. + G.G.I, + G.G.I 

V 1 J 1- 4 1 J 2. . l j 3. , i J 


D,ij S 


ref 


bJ 


i>J 


i.J 


4. 

i.J 


+ G.G.I + G. 
l 1 5 . . i 

i.J i >J' 


i.G.i, } 

1 J 6 < J 


( 6 ) 


where the G., G., G., etc., are linear functions of the v.'s (ref. 4). 

X X X ’ J 

The 1^ through 1^ are again determined by wake geometry, 

i.J i> j 

being integrals of combinations of the A, through A. times 4. n , for 

1 . . 4 . . i 

il ij 

n = 0, 1, 2, as given in reference 4. 
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The Munk optimization procedure uses (ref. 6) 


w 


— — - = w = constant 

COS<f> o 


( 7 ) 


which for the assumed wake model is written as 


COS' 


4 = ^£J( !j± ^)( a i.. + a 2 ..) 

O 3=1 \ ' 13 13 

L )(\ 4 + \ 1 )} 




13 13 


( 8 ) 


This yields N equations for N + 1 unknowns, the N y^ values and w q . 
The system is completed by specifying a value, since (ref. 4) 

I N / /Y-. n 2y..v N r (s.K) w 

{i + J-L)) * £ (“>M j s j - s -L-)J f 


C L S 


(9) 


The direct optimization procedure extremizes (minimizes) the function 
N 




( 10 ) 


where A is a Lagrange multiplier and C . is an analytical expression 

L,3 

for the derivative of equation (9) above with respect to (y./U). Similarly, 

expressions for derivatives of C D with respect to (y./U) have been 

developed, as reported in reference 4. The current method, which yields 
identical results, is to write the induced drag in matrix form as 


<*- U) T a U ) 


(ii) 


The optimal y /U values are then determined by N+l equations given by 
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= 0, i = 1, . . 


• 9 


N 


( 12 ) 



+ 



+ A C 


L, 1 


and 


N 


E 

j=i 



= 0 


(13) 


The A matrix, as given in reference 8, is in terms of the 1^ through 

i.j 

Ig from equation (6) . More details of the theory may be found in 

i, j 

references 4 and 8. 


SAMPLE RESULTS 


Convergence studies for the present method have been given in reference 
4. In general, the advanced panel method has been shown to be on the order 
of four times more accurate than a discrete vortex wake model having the 
same number of singularity unknowns, both for isolated planar wings and some 
limited, isolated nonplanar examples from references 12 and 13. In reference 
8, results are given for multiple planform configurations from references 
14 and 15. In these previous studies, the current wake model yielded induced 
drag values within 1 percent of the exact results for from 25 to 50 unknown 
wake strength values. The two optimization methods yield essentially identical 
results, except for y values near a wing tip. Cosine wake spacing greatly 
improves accuracy for a fixed number of unknowns. The reader is referred to 
references 4 and 8 for details of these studies. 

In this section some additional solutions will be presented to illustrate 
the utility of the present theory. Figure 2 displays results for the present 
theory compared with that of references 14 and 15 for the planform sketched in 
the figure. The inboard 50 percent of the wing is flat, with a constant 30- 
degree dihedral outboard of the flat portion. In addition, a vertical fence of 
variable height £ is positioned at the dihedral break span location. This is 
termed configuration 5 in reference 14. Shown in figure 2 are values of the 
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I 


induced drag efficiency parameter, k, defined as the ratio of the induced drag 
for a planar wing of equal span divided by the for the nonplanar configu- 

ration. The two theories agree favorably. For the present method, the fence has 
been oriented at <f> = 89.7 degrees to avoid numerical difficulties, as mentioned 
in reference 4. Figure 3 presents induced drag results for vee wings, compared 
with an exact solution from reference 16. In figure 4, similar results are 
presented for a diamond wake shape, again compared with exact results from 
reference 16. This last wake shape is of interest for the joined wing concept 
of reference 3. The present theory essentially duplicates these exact results. 

In figure 5 the bound circulation distributions from the present theory are 
shown for vee and diamond wings having cj) = 30°. For minimum drag, both wings 
of the diamond shaped wake carry the same lift and have the same T distribu- 
tions. The vee wing has relatively a smaller fraction of its lift developed 
inboard than does the diamond wing having the same dihedral . 

The results described above illustrate the capabilities of the computer 
program to accurately duplicate known exact solutions. As examples of more 
complicated wake geometries, figures 6 and 7 display the computed induced 
drag efficiency factors for a series of diamond wings (ref. 16), with the 
addition of end plates and winglets, respectively. These configurations can 
improve the induced drag efficiency factors for the concept of reference 3. 


CONCLUSIONS 


An advanced panel Trefftz plane wake model has been developed which 
allows accurate computation of the induced drag, bound circulation 
distribution, and wake vortex strength for nonplanar multiple planform 
configurations. The computer program which has been written to implement 
this theoretical method has been documented herein in the appendixes, 
including a listing of the code and user input instructions. A brief 
outline of the theory and some sample results have been given. These 
results reproduce accepted exact solutions for vee and diamond wings. 
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APPENDIX A 


DESCRIPTION OF COMPUTER PROGRAM AND LIMITATIONS 


This Appendix briefly describes the organization of the computer program 
written to implement the theory outlined earlier in this report, which has 
been described in some detail in references 4 and 8. Some limitations of 
this computer program are also discussed. 

This program has been written in FORTRAN IV and is currently operational 
on a Cyber 173 computer at NASA/LaRC. This computer uses approximately 15 
decimal digits in all computations. Some modifications to the code will be 
necessary if it is to be used on a computer system which uses a signifi- 
cantly different number of decimal digits. For example, the tolerances 
in subroutines SNTAN and LOGS may have to be varied. Further, double pre- 
cision arithmetic will be required for all calculations for machines using 
eight significant figures. This would entail an IMPLICIT DOUBLE PRECISION 
(A-H, 0-Z) statement in the main program and all subroutines, as well as 
use of double precision on all special functions: DCOS, DSIN, DLOG, DATAN, 

DATAN2, DSQRT, DABS, DMIN1, and DMAX1. Further, some of the variable 
names may need to be changed to be consistent with the implicit double pre- 
cision statements. 

The computer program consists of a main program, DNWASH, which reads 
the input data, performs the initial wake geometry computations, sets up 
the optimal induced drag matrix, and computes the final induced drag and 
normal wash and bound circulation distributions. This program calls nine 
subroutines: CCAL and CONCAL, which compute wake geometry constants; SNTAN 

and LOGS, which compute integrals appearing in the expression for C^, as 

detailed in the appendix to reference 4; WCAL, which computes an element 
of the optimization matrix for the Munk's criterion procedure; DRACAL and 
OPTCAL, which compute elements of the direct optimization matrix; GAMCAL, 
which computes the bound circulation terms; and SIMEQ, a linear equation 
solver. A listing of the complete computer program is given in Appendix E 
of this report, and an example input and output are given in Appendix D for 
one of the configurations discussed earlier in this report. 
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The known limitations of the computer program are now briefly described. 
First, the user-specified local dihedral angles may not anywhere equal 90 
degrees. As discussed in reference 4, this value of <j> may be approached 
(<f> = 89.5°) to approximate the wake geometry for a vertical end plate or 
pylon. (Further examples of wake geometries with nearly vertical surfaces 
have been given in figures 2 and 6.) Second, the total number of wake 
panels for all planforms is currently limited to a maximum of 50. Based 
on results from reference 4 for isolated planar and nonplanar planforms, 
this should provide an induced drag solution accuracy comparable to that 
obtainable from 200 to 250 discrete vortex unknowns; that is, better than one 
percent accuracy. Third, the maximum number of individual planforms possible 
is currently 10, while the maximum number for which runs have been attempted 
to date is only 3. Fourth, the code currently does both a Munk criterion 
optimization and a direct optimization only for a single nonplanar or 
planar planform. Solutions for multiple interacting surfaces are computed 
using only the Munk criterion solution. (It is to be noted that the design 
code described in reference 8 does have the multiple planform capability 
using the direct optimization technique.) Fifth, based upon previous 
experience (ref. 4), it is recommended that a cosine spacing of the wake 
panels on all planforms be used. 

Next, for configurations with multiple planforms, either the wakes 
must not cross one another, or any such wake crossings must occur at the 
edges of wake panels. This can be accomplished by specifying the wake 
crossing point as a common wake breakpoint on both planforms. (See 
Appendix B for a description of preparation of an input deck and definitions 
of the input data.) If wake crossings occur in the midrange of any wake 
panel, a message "80 ENTERED" is printed on the output file. For such 
cases, a midrange singularity occurs in the inverse tangent integrals 
evaluated in SNTAN. The code attempts to fit a pair of quadratics, one 
on either side of the singularity, to the inverse tangent portion of these 
integrands. The accuracy of this procedure is unknown, and any results so 
obtained are likewise of unknown accuracy. Further, for wake shapes comprised 
of continuously varying curved surfaces, it is possible for this problem of 
a midrange singularity to occur for multiple planforms even when the wakes do 
not themselves cross. Instead, all that is required to cause the apparent 
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singularity is for the projection of the plane containing one wake panel to 
intersect another wake panel away from that panel's edges. Again, this 
problem can be avoided by defining such points to be wake breakpoints on the 
second planform. It is believed that this apparent singularity is only due 
to the way in which the computer program is structured, where for example the 
above-mentioned inverse tangent integrals always occur in pairs, but each 
integral is evaluated individually. This has been alluded to in reference 4 
(p. 16), where it is remarked that integrals of the form 

s 

/ 

-s 


4 


R+T4 


tan 


-1 I c+2B4 


R+T-4 


d-4 


for 0 < n < 4, become infinite for R = T = 0. However, since what actually 
must be evaluated is an integral of the form 


/ 


R+T4 


-1 

tan 


C+2B4 

-1 

d+2B4 "| 

|R+T4 | _ 

- tan 

|R+T4 | J 


dA 


these two integrals are equal to the sum of the finite parts of the individual 
integrals, which have the form 


s x n 

-/ c+2Bi 


However, there is currently no logic in the code to automatically replace the 
original integrand with the simpler, finite part in the viscinity of a wake 
crossing point. 

Finally, the optimum wake vortex sheet strengths and bound circulations 
for a single planform for the Munk criterion solution differ slightly from 
those computed by the direct optimization solution technique. Usually these 
differences are confined to the tip region of a planform. This is believed 
to be due to the inaccuracy of the piecewise linearly varying functional form 
of the wake vortex sheet strength in the vicinity of the tip, where the 
actual wake sheet strength should be infinite. Comparisons between the two 
solution techniques for a planar isolated wing have been given in reference 4. 
These slight differences in the y and T distributions for the two 
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solution techniques do not appreciably affect the computed induced drag 
efficiency factors, but do lead to inaccuracies in the normal wash computations, 
especially for a nearly vertical surface, near a wing tip. Use of cosine wake 
panel spacing, as recommended above, will minimize this problem. 
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APPENDIX B 


INPUT DATA PREPARATION 


In this appendix the information necessary to prepare an input deck to 
use the computer program listed in Appendix E is given. A sample input deck, 
as well as the resultant output, are given in Appendix D. 

The first four input cards specify control integers and integers which 
define the number of lifting surfaces and distribution of wake vortex panels. 
These cards are all in a 515 format. The specific information needed on each 
card is as follows: 

Card Number Variable Columns Description 

1 NLLINE (1-5) Total number of lifting surfaces 

for current configuration; NLLINE 

< 10 . 

specified above, cards 2 through 5 must 


A control integer which determines 
the wake panel spacing on the 
current lifting surface, as follows: 

A. ICNTRL = 3; general input. User must 
specify NTOT (card 2, below), followed 
by values of (YHH(I) , ZHH(I), PPP(I), 
1=1, NTOT), the wake panel corner 
points and dihedral angles for the 
current surface. The YHH, ZHH, PPP 
values are specified in a 6F10.0 
format. Cards 3 to 5 described 
below are not needed when ICNTRL = 3. 

B. ICNTRL = 6; circular arc wake. User 
must input NTOT (card 2, below), 
followed by one card giving the values 
of BET, THET, and BOT in 6F10.0 format. 
BET equals the ratio of the maximum 
vertical extent of the surface to the 
semispan. A value of BET = 0. corre- 
sponds to a flat surface, while BET = 
1.0 corresponds to a semicircular 

arc. (See reference 4 for results for 
this type of surface.) THET is equal 
to the value, in degrees, for the angle 
subtended by one-half of the circular 
arc wake, while BOT equals the desired 


For each of the NLLINE surfaces 
be specified, as follows: 

2 ICNTRL (1-5) 
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Card Number 


Variable 


Columns 


Description 


wake semispan. Cards 3 to 5 described 
below are not needed when ICNTRL = 6. 

C. ICNTRL = 7; equal wake panel spacing. 

The computer program will automatically 
panel each flat portion of the wake 
between adjacent wake breakpoints 

(YY, Z) (described below; card 5), with 
equally sized wake panels. Wake panel 
size may vary on different flat portions 
of the wake. When ICNTRL = 7, card 2 
must be followed by cards 3 to 5 
described below. 

D. ICNTRL =8; cosine wake panel spacing. 
The computer program will automatically 
panel each flat portion of the wake 
between adjacent wake breakpoints 

(YY, Z) , with cosine-spaced wake panels. 
If the flat wake portion ends at the 
configuration centerline or lifting 
surface's junction with another surface, 
a quarter-circle distribution is used; 
otherwise a semicircle distribution is 
generated. This is generally the 
recommended value of ICNTRL for maxi- 
mum accuracy (see Appendix A) . When 
ICNTRL = 8, card 2 must be followed 
by cards 3 to 5 described below. 

2 NTOT (6-10) Total number of wake vortex panels on 

current aerodynamic surface. Note 
that NTOT < 50. Further, the grand 
total of the sum of all NTOT values 
for all NLLINE surfaces must not 
exceed 50. 

For values of ICNTL = 7 or 8, cards 3 to 5 must be specified as 
follows: 

3 NBRK (1-5) Total number of wake breakpoints for 

the current surface. Note that NBRK 
equals the number of changes in wing 
dihedral, plus two; or the number of 
flat portions of the surface, plus 
one. 

User-specified numbers of wake panels 
on each of the (NBRK-1) flat portions 
of the current aerodynamic surface, 
beginning at the root. 


4 LSEG(I) , 1=1, (1-25) 

.... (NBRK-1) 
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The following card is in a 6F10.0 format: 


Card Number Variable Columns 

5 YY(I), Z(I), (1-60) 

PP(I), I = 

1, .... NBRK 


After all NLLINE sets of geometry 


are needed: 

6 CLDES (1-10) 

6 SREF (11-20) 

7 - (5) 

7 - (10) 


Description 

Values of the Y, Z (in appropriate 

units), and (in degrees) for break- 

points of the current aerodynamic 
surface, beginning at the root. Note 
that the left half of the assumed- 
symmetrical planform is input, so that 
Y becomes negative going root-to-tip, 
while Z is negative up (see fig. 1). 
Note also that the PP(I) value is the 
dihedral value, in degrees, inboard 
of breakpoint I; the root value of PP 
is therefore not needed. 

data have been input, the following cards 


Desired lift coefficient, in F10.0 
format. 

Total configuration reference area 
(in appropriate units), in F10.0 
format . 

zero (0) 

zero (0) 


Card 7 signifies the end of input data for one configuration. Further 
configurations may follow card 7 beginning again with card 1. At the end 
of all configuration data for any one run, a final blank card must be included 
to signify the end of that run: 

8 - (5) zero (0) 

Note that when ICNTL equals 7 or 8 the breakpoint data specified includes 
YY, Z, and PP, which in effect overspecifies the wake geometry. This has 
not proved to be a problem, except that for more complicated configurations 
the entire F10.0 data field for each YY, Z, or PP value should contain 
significant figures for optimum accuracy of the input geometry. 
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APPENDIX C 


OUTPUT DATA DESCRIPTION 

The computer program prints out information of two general types: first, 

geometry data, both as input data and the calculated wake panel geometry, are 
printed. This information is followed by the minimum drag solution informa- 
tion, which includes the wake vortex sheet strengths, optimum bound 
circulation, induced drag coefficient, and induced drag efficiency factor. 

For a single planform, this solution information is printed for the Munk's 
criterion solution, followed by the same output for the direct optimization 
technique solution, while for configurations with more than one planform, 
only the Munk's criterion solution is computed and printed. In this 
Appendix, the output information for a configuration is described in the 
order in which it printed. 


Geometry Data 


For each lifting surface of the configuration the values of the wake 
breakpoints, from root to tip (see Appendix B) , are listed. This is 
followed by the peripheral length of that surface. Next, the individual 
wake vortex panel corner points, dihedral angles (in radians), and panel 
semiwidths are listed for the entire configuration. Finally, reference 
quantities for the configuration are listed. In detail, the geometry data 
listed is as follows: 


(YY (I) , 1=1, . . . , NBRK) 

(Z(I) , 1=1, .... NBRK) 
(PP(I) , 1=1, . . . , NBRK) 


Y coordinates of the breakpoints of each 
lifting surface, as input by the user, 
ordered from root to tip. (See figure 1 
for positive coordinate directions.) 

Z coordinates of the breakpoints of each 
lifting surface, as input by the user, 
ordered from root to tip. 

Dihedral angles, <f>, just inboard of break- 
point I, in radians. 


DTOT 


Total peripheral length of each lifting 
surface . 
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The above data is followed by the wake panel corner points and semiwidths 
actually used in the Trefftz plane calculations. First, the program lists 
whether equal or cosine spaced wake panels have been generated. This is 
followed by a table containing the following information: 


i; i = i. 


(NTOTT +1) 


YHH(I) 

ZHH(I) 

PPP(I) 


SNN(I) 


Individual wake panel number, numbered 
from the tip of the first lifting surface 
to the root, followed by the tip-to-root 
numbering of wake panels on successive 
surfaces. NTOTT is the total number of 
wake panels for the entire configuration; 
NTOTT < 50. 

Y coordinate of outboard, or tipmost, 
corner of wake panel I. 

Z coordinate of outboard corner of wake 
panel I. 

Dihedral angle, in radians, of wake panel 

I. 

Semiwidth of wake panel I. 


Finally, the following reference quantities are listed: 


T0L2 


CLDES 

SREF 

BSAVE 

ARAT 


The tolerance utilized in subroutine LOGS . 
Generally, the value of the tolerance 
utilized in SNTAN, TOL, will have the 
same value, unless changed by the user. 

The value of T0L2 should be small compared 
to the smallest value of SNN. 

The desired lift coefficient value. 

The configuration reference area. 

The configuration reference span, taken 
as twice the maximum absolute value of 
YHH. 

The configuration aspect ratio, defined 
as (BSAVE) 2 /SREF. 


Solution Data 


The output data for the minimum drag solution for the Munk's criterion 
solution consists, first, of a table of the following: 
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I Individual wake panel number, numbered 

tip-to-root, as described above, for each 
lifting surface of the configuration. 

BGAM(I) Bound circulation value, r/U, at the out- 

board, or tipmost, corner of a wake panel 
I for minimum induced drag at a specified 
lift coefficient. 

CDRAG(I) Nondimensional bound circulation value 

for minimum induced drag, Y/Y, at the 
outboard corner of wake panel I. Note 
that for a wake consisting of a portion 
of a circular arc, CDRAG values are non- 
dimensionalized by the F value at the 
root of the planform. 

AOPT(I, NTOTT + 1) Wake vortex sheet strengths for minimum 

induced drag, y/U, at the outboard 

corner of wake panel I. 

GAM(I) Nondimensional wake vortex sheet strength 

for minimum induced drag, y/y, at the 
outboard corner of wake panel 1. 

ETA Nondimensional spanwise coordinate of 

outboard corner point of wake panel I, 
at which the above values are computed. 

It is after this information, during the computation of the induced drag, 
that it is possible that a message "80 ENTERED" may be printed, to indicate 
a problem with wakes crossing one another, as discussed in Appendix A. This 
is followed by the induced drag coefficient, induced drag efficiency factor, 
and computed normal wash velocities as follows: 


CD 


DIDEAL 


WDBU 


DEFF 


Induced drag coefficient computed using 
optimum value vortex strengths, for Hunk's 
criterion solution. 

Induced drag coefficient for a planar wing 
having the same projected span as the 
current configuration, evaluated at the 
same lift coefficient value. 

The ratio of the constant, w Q , appearing 
in the general statement of Munk's criterion 
[eq. (7)], divided by U. 

Induced drag efficiency factor, k, for 
the configuration, defined as the ratio 
of the induced drag for the planar wing 
divided by the computed induced drag for 
the configuration. 
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I 


Wake panel number. 


WDOWN Computed induced normal velocity at the 

midpoint of wake panel I. 

WOP Induced normal velocity divided by the 

cosine of the dihedral angle, evaluated 
at the midpoint of wake panel I. 

CDAPP An approximate value of induced drag 

coefficient, evaluated by assuming T 
and y are constant on each wake panel. 

For single planforms, all of the above output, with the exception of the 
initial geometry data, is repeated for a second solution achieved using the 
direct optimization procedure for the same configuration. A sample output, 
as well as the input data deck, appear in Appendix D. 
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APPENDIX D 


EXAMPLE OF INPUT AND OUTPUT DATA 


Sample input data and output data are presented for one of the configu- 
rations of figure 6 of this report, were d/h = 1.0, h/b = 0.355. Input data 
and a sketch of the input wake shape appear on page 22 , while the output 
data begins on page 23. 
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UJ cd r\) 


Input Data and Sketch of Wake for 
Diamond Wing with End Plates; 
d/h = 1.0, h/b = 0.355 


15 

1 0 5 

*0.1775 0.0 -0.4-98451 0.0 -19.600981 

-.5 -.17750 89.5 

8 15 

3 

10 5 

°* 0.1775 0.0 -0.498451 0.0 19.600981 

-.5 .17750 -89.5 

0.5 1.0 

0 0 

0 
0 
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Output Data for Diamond Wing with 
End Plates; d/h = 1.0, h/b = 0.355 


GENERAL INPUT GEOMETRY 


0.00000 

-.49845 

-.50000 

TOTAL PLANFORM 


-.17750 0.00000 

0.00000 -19.60098 

• -.17750 89.50000 

PERIPHERAL LENGTH= 


.70662 


0.00000 

-.49845 

-.50000 

TOTAL PLANFORM 


.17750 0.00000 

0.00000 19.60098 

.17750 -89.50000 

PERIPHERAL LENGTH= 


.70662 


COSINE SEGMENT SPACING 


SEGMT NO 

Y 

Z 

PHI 

1 

-.499926 

-.169025 

1.562070 

2 

-.499658 

-.138363 

1.562070 

3 

-.499225 

-.088750 

1.562070 

4 

-.498793 

-.039137 

1.562070 

5 

-.498525 

-.008475 

1.562070 

6 

-.495383 

-.001093 

-.342102 

7 

-.483185 

-.005436 

-.342102 

8 

-.459089 

-.014017 

-.342102 

9 

-.423689 

-.026623 

-.342102 

1 0 

-.377857 

- .042944 

-.342102 

1 1 

-.322720 

-.062578 

-.342102 

12 

-.259637 

-.085042 

-.342102 

13 

-.190161 

-.109783 

-.3421 02 

1 4 

-.116002 

-.136191 

-.342102 

15 

-.038987 

-.163616 

-.342102 

16 

-.499926 

.169025 

-1.562070 

17 

-.499658 

.138363 

-1.562070 

ie 

-.499225 

.088750 

-1.562070 

19 

-.498793 

.039137 

-1.562070 

20 

-.498525 

.008475 

-1.562070 

21 

-.495383 

.001093 

.342102 

22 

-.483185 

.005436 

.342102 

23 

-.459089 

.014017 

.342102 

24 

-.423689 

.026623 

.342102 

25 

-.377857 

.042944 

.342102 

26 

-.322720 

.062578 

.342102 

27 

-.259637 

.085042 

.342102 

28 

-.190161 

. 1 09783 

.342102 

29 

-.116002 

.136191 

.342102 



Output Data for Diamond Wing with End Plates; 
d/h = 1.0, h/b = 0.355 (continued) 


30 -.038987 .163616 



SNN(I) 

1 

.00848 

2 

.02219 

3 

.02743 

4 

.02219 

5 

.00848 

6 

.00326 

7 

.00969 

8 

.01589 

9 

.02169 

10 

. 02696 

11 

.03157 

12 

.03540 

13 

.03835 

14 

.04037 

15 

.04139 

16 

.00848 

17 

.02219 

18 

. 02743 

19 

.02219 

20 

.00848 

21 

.00326 

22 

.00969 

23 

.01589 

24 

.02169 

25 

.02696 

26 

.03157 

27 

.03540 

28 

.03835 

29 

.04037 

30 

.04139 


TOL IN SNTAN = . 26057E-05 


DESIGN LIFT COEF. = .50000 

WING REFERENCE AREA = 1.00000 

REF WING SPAN = 1.00000 

ASPECT RATIO = 1.00000 


.342102 
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Output Data for Diamond Wing with End Plates; 
d/h = 1.0, h/b = 0.355 (continued) 

OPTIMUM LOADING USING MIJNKS CRITERION 


SEGMT 

BOUND CIRC 

bgam/ave 

SNED STRH , 

GAM/AVE 

ETA 

1 

0. 

0 . 

■ 3 09 0 7E + 0 1 

.31321E+02 

. 10000E+01 

2 

• 33090E-0 1 

•29612E+00 

, 8 1 359E + 0 0 

.82448E+0 1 

. 99970E + 0 0 

3 

, 62488E-0 1 

•55920E+00 

.51 136E+00 

.518P0E+01 

.99893E+00 

4 

.8 7560E-0 1 

.78357E+00 

•40281E+00 

.40820E+01 

.99797E+00 

5 

• 1 0456E + 0 0 

. 9 3567E + 0 0 

. 36 32 0E + 0 0 

• 368 06E + 0 1 

.99720E+00 

6 

. 1 0955E+00 

,98 0 35E + 00 

• 22603E+00 

. 22905E+ 0 1 

. 9969 OE + 0 0 

7 

• 1 1 04 OE + 0 0 

. 98795E + 0 0 

.34701E-01 

. 35 1 65E + 0 0 

.98463E+00 

8 

.11089E+00 

.99232E+00 

. 15637E-01 

. 15847E+00 

. 948 1 1 E * 0 0 

9 

. 1 1 156E+00 

,99831 E+00 

• 26488E-0 1 

. 2684 2E + 0 0 

.88825E+00 

1 0 

• 1 1 3 22E + 0 0 

. 10132E+01 

.50322E-01 

.50996E+00 

.80651E+00 

1 1 

. 1 1658E+00 

, 10432E+01 

.74120E-01 

.751 13E+00 

. 70492E+00 

12 

. 12160E+00 

. 10862E+01 

. 84958E-0 1 

.86096E+00 

.58596E+00 

13 

. 12748E+00 

. 1 1408E+01 

.S1092E-01 

. 82 1 78E + 0 0 

. 45 25 8E + 0 0 

14 

. 1331 4E + 00 

. 1 1915E+01 

. 66558E-0 1 

.67449E+00 

. 30806E+00 

15 

. 13762E+00 

. 12315E+01 

• 444 1 OE-O 1 

.45004E+00 

. 15595E+00 

16 

0. 

0. 

.30907E+01 

• 3 1 3? 1 E + 02 

.10000E+01 

17 

•33090E-01 

. 296 1 2E + 0 0 

. 8 1 359E ♦ 0 0 

. 8244 fiE + 0 1 

. 9997 OE ♦ 0 0 

18 

. 62488E-0 1 

.55920E+00 

.51 136E+00 

.51820E+01 

. 99893E+00 

19 

. 8 756 0E- 0 1 

. 78357E+00 

•40281E+00 

.40820E+01 

.99797E+00 

20 

. 10456E+00 

. 93567E + 0 0 

. 36320E+0Q 

. 368 06E + 0 1 

. 99720E+ 00 

21 

. 10955E+00 

. 980 35E + 0 0 

•22603E+00 

. 22905E + 0 1 

.996R0E+00 

22 

. 1 1 040E+00 

, 98795E + 00 

.34701E-01 

. 35 1 65E + 0 0 

.98463E+00 

23 

. 1 1089E+00 

.99232E+00 

. 15637E-0 1 

. 15847E + 00 

.9481 1E+00 

24 

• 1 1 1 56E + 0 0 

.998 3 1 E + 0 0 

.26488E-01 

. 26842E ♦ 00 

.88825E+00 

25 

. 1 1322E+00 

. 10132E+01 

. 50322E-0 1 

.50996E+00 

.80651E+00 

26 

. 1 1658E+00 

. 10432E+01 

•74120E-01 

.751 13E+00 

.70492E+00 

27 

. 1216 OE *00 

. 1 0882E + 01 

. 84958E-0 1 

.86096E+00 

.58596E+00 

28 

, 12748E+00 

. 1 1408E+01 

.81 092E-0 1 

.82178E+00 

, 45258E + 0 0 

29 

. 13314E+00 

.1 1915E+01 

. 66558E-0 1 

. 674 49E + 0 0 

.30806E+00 

30 

. 1 376 2E ♦ 0 0 

. 12315E+01 

.4441 0E-01 

.45 0 0 4E + 00 

. 15595E + 00 

31 

. 13946E+00 

. 124P.0E + 01 

0. 

0. 

0. 


200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTEREO 
200 ENTERFD 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
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Output Data for Diamond Wing with End Plates; 
d/h = 1.0, h/b = 0.355 (continued) 


200 ENTERED 

200 ENTERED 
200 ENTERED 
200 ENTERED 

POO entered 

2O0 ENTERED 
200. ENTERED 
200 ENTERED 
20 0 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 
200 ENTERED 


CD CALCULATED USING SUB DRACAL AND OPTIN LOADS USING MUNK CPIT = 

.47846E-01 

CD FOR FLAT WING = .79577E-01 

RATIO OR ZERO DIHEDRAL DOwNwSH/U = .38317E+00 


INDUCED DRAG EFFICIENCY 


FOR WINGS OF EQUAL SPAN 


. 166319E+01 


1 

DOWNWASH 

I 

.83594E-0 

2 

, 83594E-0 

3 

.83594E-0 

4 

. A3B94E-0 

5 

.B3594E-0 

6 

, 9024?E-0 

7 

. Q 0P4?E-0 

8 

.90242F-0 

9 

.90P42F -0 

10 

. 90 24 ?E-0 

1 1 

.90242E-0 

12 

. 9024?E-0 

13 

, 90242E-0 

14 

.90242E-0 


W/COS {PHI ) 

3 . 9 5793E - 0 1 
3 .95793E-01 
3 . 95793E-0 1 
3 . 95793E-0 1 
3 .95793E-01 
1 • 95793E- 0 1 
1 • 95793E-0 1 
1 . 95793E-0 1 
1 .95793E-01 
1 . 95793E- 0 1 
1 .95793E-01 
1 .95793E-0 1 
1 .95793E-01 
1 .95793E-01 
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I 


Output Data for Diamond Wing with End Plates; 
d/h = 1.0, h/b = 0.355 (concluded) 


15 

• 90 242E-0 1 

.95793E-01 

16 

.83594E-03 

.95793E-01 

17 

, 83594E-03 

, 95793E-0 1 

IP 

.83594E-03 

.95793E-01 

19 

.P3594E-03 

• 95793E-0 1 

20 

.83594E-03 

.95793E-01 

21 

.90242E-01 

.95793E-01 

22 

■ 902426-0 1 

.95793E-01 

23 

,90 24 ?E- 0 1 

.95793E-0 1 

24 

.90242E-01 

.95793E-01 

25 

.90242E-0 1 

. 95793E-0 1 

26 

, 90242E-0 1 

.95793E-01 

27 

, 90242E-0 1 

.95793E-01 

28 

.90242E-0 1 

, 957 93E-0 1 

29 

.90242E-01 

. 95793E-0 1 

30 

. 9Q242E-0 1 

, 95793E-0 1 


APPPOX CD USING SOLVED BOUND 

CIRCULATIONS AND WASHES AT SEG MIDPOINTS = .47909E-01 
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APPENDIX E 


COMPUTER PROGRAM LISTING 

This program has been written in FORTRAN IV language for the CDC series 
6000 computer system with N0S1.3 operating system. Minor modifications may 
be necessary to achieve successful execution on other computers, as discussed 
in Appendix A. The following table is an index to the computer program 
listing : 


Name 


Letter Designation 

Page 

PROGRAM DNWASH 

A 

29 

SUBROUTINE 

GAMCAL 

B 

47 

SUBROUTINE 

WCAL 

C 

48 

SUBROUTINE 

CCAL 

D 

50 

SUBROUTINE 

CONCAL 

E 

51 

SUBROUTINE 

SNTAN 

F 

52 

SUBROUTINE 

LOGS 

G 

60 

SUBROUTINE 

DRACAL 

H 

63 

SUBROUTINE 

SIMEQ 

I 

65 


The permanent file name of this program at NASA/Langley Research 
Center is DRG, stored under user number 496125E. 
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PROGRAM PNW4SH 


PAGE 


l 


1 


S 


10 


15 


?0 


?5 


30 


35 


40 


C 

C 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 


PROGRAM PNW4SH ( INPUT .Oil I PI 11 , T APES « T APE 6 = OUTPtlT ) 

0 1 PENS I ON GARISH. HGAM(Sl). VHEUSl), 7HH(51), PPPI51) 

D I ME MS 1 ON PNNSH(Sl) 

DIMENSION AOPT (52,5?) 

DIMENSION AIM (6) 

DIMENSION COR AG (51) 

DIMENSION ASlM(Sl.Sl). PSIMI51.1), IPIV0T(51) 

DIMENSION YY(IO), 7(10). PP(10). LSEG(IO), DTHETE (10) 

DIMENSION PEPIF(IO), H ( 1 0 1 . GSUM(IO) 

DIMENSION CLP ( 53 ) » SEO(IO). 010(10) 

COMMON /TFLL/ TOI.TOL7.IOL3 
COMMON /SEG/ SNN(Sl) 

COMMON /DIPOPT/ T 1 (53.531.12(53.531.13(53.53) .TA(53.53).T5(53,53) 


1,16(53.531 

COMMON /f EN/ NSPKlOl.NLl 1NE 
PI = 4.«ATAN(1.) 

DTP = PI/1H0. 

PROGRAM ViPITTFN TO IMPLEMENT TRtEFT7 PLANE INDUCED DRAG 
OPTIMIZATION THEORY DESCRIBED IN NASA CR-3154, JUNE 1979. 


PROGRAM WRITTEN BY DR. JOHN M, KUHLMAN, DEPT MECHANICAL 
ENGINEERING AND MECHANICS, OLD DOMINION UNIVERSITY, NOP- 
EOLK, VA 23S0H. UNDER NASA LANGLEY GRANT NSG-1357. 

THEORY ASSURES A TWO DIMENSIONAL ADVANCED PANFL MODEL OF 
THE UNO I STORIED , INTERACTING WAKES OF MULTIPLE LIFTING 
SURFACES. 


WAKE VORTEX SHFFT STRENGTHS ARE ASSUMED TO VARY IN A 
PIFCEWISE LINEAR FASHION. ANALYTICAL EXPRESSIONS FOP 
INDUCFD NORMAI VFLOCITY, BOUND CIRCULATION, INDUCED DRAG, 

AND LIFT ARE DEVELOPED IN CR-3154 IN TERMS OF THE ASSUMED 
WAKE MODEL. 

THESE EXPRESSIONS ARE USED TO OBTAIN MINIMUM DRAG WAKF VORTEX 
SHFFT STRENGTHS. BOUND CIRCULATION DISTRIBUTIONS, AND CD 
VALUFS for minimum DRAG »T A GIVEN LIFT, FOP NONPLANAR MULTIPLE 
INTERACIING LIFTING SURFACES, USING BOTH MUNK • S CRITERION AND 
A 0 1 RFC! OPT I M I 7 A T I ON PROCEDURE. 


ICTRL = 3 FOR GENERAL INPUT Of NON-Pl.ANAR WING CASES 


A 1 
A 2 
A 3 
A 4 
A 5 
A 6 
A 7 
A 8 
A 9 
A 10 
A 11 
A 12 
A 13 
A 14 
A 15 
A 16 
A 17 
A 18 
A 19 
A 20 
A 21 
A 22 
A 23 
A 24 
A 25 
A 26 
A 27 
A 28 
A 29 
A 30 
A ‘31 
A 32 
A 33 
A 34 
A 35 
A 36 
A 37 
A 38 
A 39 
A 40 
A 41 
A 42 


w 

>£> 


u > 

o 


PROGRAM DNWASH 


PAGE 


3 



C 


TMFN MUST INPUT THIPIE5 OF YHH.7HH.PPP FOP EACH 

CORNER POINT 


A 

43 


C 


OF EACH SEGMENT, STARTING AT WING 

TIP 





A 

4* 

AG 

c 


I C TPL = 6 FOR CIRCULAR ARC RING, CONE, 

TR 1 39 




A 

45 


c 


ICTRL=7 FOR WING OF ARBITRARY PHI 

AND 

ETA-EOUAL 

SPACING 


A 

46 


c 


lCTPL=fl FOH WING OF arbitrary PHI 

AND 

f TA-COSINF 

SPACING 


A 

47 


c 








A 

48 


c 


HTOT=NUMHFR OF SEGMFNTS ON SfNISPAN, 

BOTH FOR FOUAL 

AMO COSINE 

SPA 

A 

49 

GO 

c 








A 

50 


c 


REAR AN INTEGFP VARIABLE Nl.LINE = 

BENT 

LIFTING LINES 

NEEDED TO 

MAKE 

A 

5) 


c 


NTOTT = NUMbFR OF SFGMENTS ON NLLINF. 





A 

53 


c 


NS = START NUMBER ON EACH NLLINfc 






A 

53 


f 








A 

5* 

66 


1 

RF Al) IG.ni) NLI.INE 






A 

55 




IF (niline.eo.o) go to i?q 






A 

56 




NTOTT = 0 






A 

57 




LSTART = 1 






A 

58 




DTOH = 0.0 






A 

59 

60 



WRITE (6,130) 






A 

60 




WRITE (6,136) 






A 

61 



? 

CONTINUE 






A 

63 




RF AO ( G , 1 3 1 ) ICNTHL.NTOT 






A 

63 




IF (ICNTRL.FO.O) GO TO 1 






A 

64 

66 



NS = NTOTT ♦ 1 






A 

65 




NSPT ( ES T ART ) = NS 






A 

66 




NTOTT = NTOTT.NTOT 






A 

67 




1ST OP = 1 






A 

68 




1 DRAG = 1 






A 

69 

70 

c 








A 

TO 




IF ( ICNTPL .F0.3) GO TO 13 






A 

71 




IF (ICNTPL. FO, 6) GO TO 15 






A 

73 




IF ( ICNTRL.FO.7l GO TO 30 






A 

73 




IF (ICNTPL. EG. 81 GO TO 30 






A 

74 

7G 

c 








A 

75 


c 








A 

76 


c 


CALCULATE OOWNVJASH AT SEGMENT MIDPOINTS 




A 

77 


c 








A 

78 



3 

CONTINUE 






A 

79 

80 



WRITE (6,133) 






A 

80 




00 11 1=1, NTOTT 






A 

81 




WOOWN = 0. 






A 

83 




no 10 J=l, NTOTT 






A 

83 




S = SNN(J) 






A 

BA 


PROGRAM ONWASH 


PAGE 


3 


85 CALL CCAL ( I , J,YHH, 7 HH.PPp,S*AA,U 8 ,nU,FF ,GG , EE . A J , AK , PP ♦ T T ,UU, WW) A 85 

CALL CONCAL (AA. BH, FF, GO, S, A, H,C,DfF,G,CJ,CK,CL,CM,CN, CO, CP, 1) A 86 

IF (RR.EO.O.) GO TO 4 A 8T 

P = 7.® (ATAN(C/ARS(RRH -ATAN (O/AHStRR) ) >/(ARS(RRI) A BB 

GO TO 5 A 89 

90 4 CONTINUE A 90 

P = 2./(FF-?.®S)-2./(FF*7.«S) A 91 

5 CONTINUE A 92 

IF (UU.EO.0,1 GO TO 6 A 93 

0 = 2.® (ATAN( <AJ,7.®S) /ARS(UU) I-ATAN ( (AJ-p.®S>/ABSftJU) ) I /AHSUIU) A 94 

95 GO TO 7 A 95 

6 CONTINUE A 96 

0 = 2./<AJ-?.*S)-2./(AJ,2,®S) A 97 

7 CONTINUE A .98 

Z1 = (S®S«FF»S*GG) / (S«S-FF»S*GG) A 99 

100 Z 2 = (S»S*AJ«S*AK I / (S®S-AJ»S*AK) A 100 

A 1 T J = (P*A*,S*HR®ALOG(71) ) / ( 2 , ®PI ) A 101 

A31J = (f L«P*?.«BR*S»CO» A lOG (Z1))/(2«»S*PI) A 102 

CALI CONCAL (0D,EF.AJ,AK,S,A,H,C,D,F,G,CJ t CK,CL»CM,CN,C0,CP,2> A 103 

A? I J = - (0»A, ,5«EE®AI OG(72| )/(2.®PI) A 104 

105 A4 I J = - CCL*0*2.*EE*S«C0»AL0G (Z21)/(2.*S®PI) A 105 

DO B K=1,NLLINE A 106 

KF = K ♦ 1 A 107 

KCHK = NSPT(KK)-1 A 108 

IF (J.FO.KCHK) GO TO 9 A 109 

HO B CONTINUE. A 110 

WOOWN = WOOWU« , 5® (GAN ( J* 1 1 *GAM ( J) 1 » ( A 1 1 J+ A2 1 Jl ♦ «5* (GAM | J* 1 ) -GAM ( J ) A 111 

1 1 * ( A3I J* A4 I J) A 112 

GO TO 10 A 113 

9 CONTINUE A 114 

115 WOOWN = WOOWN*. 5®GAM(J)»(A1 IJ*A2IJ-A3IJ-A4IJ) A 115 

10 CONTINUE A 116 

C A 117 

C OOVNWASH AT WI6G IS .5 HASH AT MINUS INFINITY A 118 

C A 119 

120 WOOWN = WOOWN/?, A 120 

WOP = WOOWN/COS (PPP (II) A 121 

WRITE (6,134) I, WOOWN, WOP A 122 

ONWSH(I) = WOOWN A 123 

11 CONTINUE A 124 

125 IF ( I STOP • GF , 2 ) GO TO 124 A 125 

1 ST OP = 1 ST OP ♦ 1 A 126 



w 

N> 


PROGRAM DNWASH 


PAGE 


A 


GO TO H A 127 

C A 128 

C GENERAL GEUMFTPY CAICUIATIONS A 129 

130 C A 130 

1? CONTINUE' A 131 

RF AO (“1.136) (YHH(I ) ,7HH(I) .PPPlIl ,I=NS.NTOTT| A 132 

HOT = ABS(YHH(NS) M SUN ( NS ) «C05 t PRP (NS > ) A 133 

H (l ST ART ) = ?.«HOT A 134 

135 OTOT =0. A 135 

00 13 I=1.NT0T » 136 

OTOT = OTOT + SNN ( I ) A 137 

13 COMT I Ml It A 13B 

OTOT = 2.«DTOT A 139 

140 WRITE (6.13R) OTOT A 140 

OTO(LSTART) = OTOT A 141 

OTOB = OTOFUOTOT A 142 

14 CONTINUE A 143 

GO TO 19 A 144 

145 C A 145 

C GEOMETRY FOR CIRCULAR ARC AIRFOIL A 146 

C A 147 

15 CONTINUE A 148 

MTOT1 = NT OT ♦ 1 A 149 

150 RE AO (5.136) PET.THFT.HOT A 150 

THET = THE T »PTP A 151 

R(LSTAPT) = ?.«HOT A 152 

D = HF T “HOT A 153 

R = HOT/SlN(THFT) A 154 

155 01HFTA = THFT/ (FLOAT (NTOT) ) A 155 

00 16 I=1,NT0T A 156 

YT = -R«SIN(THE T-FI.OAT ( 1-1) "OTHtTA) A 157 

2T = -R*R«C05 (THFT-FLOAT ( 1-1 ) *>DTHETA) A 158 

YT 1 = -R*S1N(TRET-FI.0AT ( l I »0THET A ) A J59 

160 2 T 1 = -RtROCOS ( THET-FLOAT ( I ) «DTHETA) A 160 

II = NS- 1*1 A 161 

PPP(II) = atari ( 7T-7T1 )/(YT-YTl> 1 A 162 

SNN(II) = .S»SORT ( (YT-YT1)*»2* (ZT-7T1) »«2) A 163 

YHH(II) = .5»(YT«YT1) A 164 

165 7HHU1) = .5* ( 7T *7T 1 ) A 165 

16 CONTINUE A 166 

OTOT = 0. A 167 

00 17 JsUS.NTOTT A 168 
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iro 


175 


ion 


185 


190 


195 


?on 


205 


?10 


DTOT = OTOI *SNN ( I ) 

17 CONTINUE 

DTOT = ?.«OTOT 
HRITt (ft, 138) DTOT 
OT 0 (LSTAPT ) = DTOT 
OTOH = OTOruOTOT 
IF (NS.NF .1 ) GO 10 18 
DTO 1 = DTOT 

18 CONTINUE 

19 CONTINUE 
GO TO 29 

C 

20 CONTINUE 
C 

C GEOMETRY CALCS FOR WAKE MADE OF STRAIGHT SEGMTS (CONCTD) 

C 

HEAD (5 •. 1371 NBRK 
NFlP = NPPK-1 

PE AO (9*137) (LSFG(I) ,I = 1,NBP) 

C NHPK FOIIALS NUMHER OF DIHEDRAL CHANGES OR JCTS WITH 
C OTHER LIFTING LINES *2 

READ (5.13ft) (YY ( I ) .7 (I ) .PP ( I ) .1=1 *NHRK) 

DO 21 1=1 ,NBRK 

21 WRITE (ft. 132) YY ( 1 ) .7 (I ) .PP( 1 ) 

DO 22 1=1 . NOP 

22 PERIF(T) = SQPT ( (7 ( 1*1 )-7 (I) ) »«2* ( Y Y ( I + 1 >-YY( I) )««2) 
DTOT = 0. 

DO 23 1=1, NBR 

23 DTOT = OTOT .PER 1 F ( I ) 

WRITE (ft. 138) OTOT 
OTO(LSTAPT) = OTOT 
DT OH = 0TO6.DT0T 

IF (NS.Nf . 1 ) GO TO 24 
DTO 1 = OTOT 

24 CONTINUE 

DO 25 1=1 ,NGRK 

25 PP ( I ) = DTR*PP ( I ) 

H(LSTART) = ?.«ARS (YY (NRRK) ) 

IF ( 1CNTRI .EO.El GO TO 30 
DO ?ft 1=1, NBR 

SEC ( I ) = 0.5°FER1F (I)/LSE6(I) 

26 CONTINUE 


A 169 
A 170 
A 171 
A 172 
A 173 
A 174 
A 175 
A 176 
A 177 
A 1 78 
A 179 
A 180 
A 181 
A 182 
A 183 
A 184 
A 185 
A 186 
A 187 
A 188 
A 189 
A 190 
A 191 
A 192 
A 193 
A 194 
A 195 
A 196 
A 197 
A 198 
A 199 
A 200 
A 201 
A 202 
A 203 
A 204 
A 205 
A 206 
A 207 
A 208 
A 209 
A 210 


UJ 

u> 



U) 

•P- 
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S = SEQ (NE)R ) 

A 

211 




SNN(NS) = S 

A 

212 




PPP(NS) = PP(NRRK) 

A 

213 




YHH (MS ) = YY (NBPK) *S»COS (PPP (NS) ( 

A 

214 

215 



7HHINS) = 7 (NtJPK) ♦S»SIN(PPP (NS) ) 

A 

2)5 




Ml = NS 

A 

216 




00 20 J=1,NBR 

A 

217 




NSEG = N1 ,LSFG (NRRK-J) -1 

A 

218 




11 = Hl.l 

A 

219 

220 



LL = NRP-J 

A 

220 




LH = LL ♦ 1 

A 

221 




00 27 1=11, USES 

A 

222 




SNN(I) = SEO(LH) 

A 

223 




YMH(l) = YHH(1-1)‘2.»SE0(LH)«C0S<PPP(M>) 

A 

224 

225 



7HH(1) = 7HHU-1 ) ♦2.*SE0(LH)«S1N(PPP(N1 ) ) 

A 

225 



27 

PPP(l) = PPP(M) 

A 

226 




Ml = Nl*t SFG(NBRK-J) 

A 

227 




IF (LL.EO.O) GO TO 20 

A 

228 




SUM (Nil = SFO(LL) 

A 

229 

230 



YHH (Ml ) = YHH(NSFG) .COS <PP (LH) ) »5EQ (LL) *COS (PPP (NSEG) >«SEQ(LH) 

A 

230 




2HH (Ml ) = 7HH(NSEG)*5IN(PP(LH) ) "SEOUL > *SIN (PPP (NSEG) )»SEQ(LH) 

A 

231 




PPP (Ml ) = PP (LH) 

A 

232 



28 

CONTINUE 

A 

233 



29 

CONTINUE 

A 

234 

235 



GO TO 42 

A 

235 


C 



A 

236 


C 


COSINE SPACING CM CULAT10NS 

A 

237 


C 



A 

23B 



30 

CONT INUE 

A 

239 

24 0 



IF 1 Y Y ( 1 ) .10,0.) GO TO 32 

A 

240 




00 31 1=1, NPR 

A 

241 



31 

OTHFTF ( I ) = PI/(FL0AT (LSFG(NRRK-I))) 

A 

242 




GO TO 35 

A 

243 



32 

00 34 1=1 ,NPP 

A 

244 

245 



IF (I.FU.NhP) GO TO 33 

A 

245 


C 


THE VARIAHLF OTHFTA IS NAMED ARRAY DTHETF(I) IN ICNTRL=« 

A 

246 




OTHFTF ( 1 ) = Pl/IFLOAT (LSEG(NBRK-l))) 

A 

247 




GO TO 34 

A 

248 



33 

CONTINUE 

A 

249 

250 



DTHFTF(I) = PI/ (2. "FLOAT (LSEG (NHHK-I) ) ) 

A 

250 



34 

CONTINUE 

A 

251 




IF (NRP.NE. . 1) GO TO 35 

A 

252 
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SNN(NSI = 0.5«PtRIF (l|#(l.-rOS(OTHETf (l( U 

A 

253 



GO TO 36 

A 

254 

255 

.75 

PH = 0.5»PFRIF (NHp) 

A 

255 



SNN(NS) = 0.5»BH»(1 .-COS (PTHETE (1) ) ) 

A 

256 


36 

PPP (N5) = PPIF.HRK) 

A 

257 



NN = 0 

A 

256 



no 40 1=1, NPP 

A 

259 

260 


N N = NNHSFGIMHPK-I ) 

A 

260 



IF (I .FO.NRP.ANn.YY U > .FO.O.) GO TO 37 

A 

261 



RH = n.5«PEPIF ( NHPK - I ) 

A 

262 



GO TO 36 

A 

263 


37 

CONTINUE 

A 

264 

265 


HH = PFRIF (NPRK-I ) 

A 

265 


3P 

CONTINUE 

A 

866 



LL = LSFG (NftpK-I ) 

A 

267 



00 39 J=l,t.L 

A 

266 



IF (I.FO.l.ANO.J.Fa.l) GO TO 39 

A 

269 

270 


MM = NN-LL«J 

A 

270 



KM = MS .NN- 1 

A 

271 



PPP(KM) = PP (NHRK . 1 - 1 ) 

A 

272 



SNN(KM) = 0,50*RH* (COS (FLOAT ( J) «OTHE TE ( 1 ) ) -COS ( FLOAT ( J- I) «DTHE TE ( I 

A 

273 



III) 

A 

274 

275 


SNM(KM) = ARS (5NN (KM) ) 

A 

275 


39 

CONTINUE 

A 

276 


AO 

CONTINUE 

A 

277 



YHHCN5) = YY (NHRK) .SNN(NS) «C0S(PPP(NS) ) 

A 

270 



7HH(NS) = 7{MIPK)«SNN(NS)»SIN<PPP<MS) ) 

A 

279 

260 


N5EG = NTOT 

A 

260 



00 41 1 1 =2 ,MSF G 

A 

281 



I = N5- 1*11 

A 

282 



YHh(I) = YF<H( 1-1 ) .SNFUl-l )»C0S(PPP(I-1 ) ) .SNN< I)»COS(PPPtI ) ) 

A 

283 



ZHHIII = 7HH(l-|)»5HN(I-l>«5IN(PPP(I-n>*SNN(I>«SIN(PPP(In 

A 

284 

285 

4] 

CONT IMIIF 

A 

285 


4? 

CONTINIIt 

A 

286 



IF (L5TART.FQ.NLLINF) GO TO 43 

A 

287 



(.START = I S T APT ♦ 1 

A 

288 



GO TO 2 

A 

289 

290 

4 3 

CONT I NUF 

A 

290 



IF UCNTPl .Fq.f ) WRITE (6.139) 

A 

291 



IF UCNTPl ,F0. 7) WRITE (6.140) 

A 

292 



IF (JCMTRL.FQ.P) WRITE (6,141) 

A 

293 



LN = MI L INF . 1 

A 

294 


LO 

Ln 



u> 

ON 
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?95 



NSPT(IN) = NTOT T + 1 

A 

?95 




WRITE (6,14?) 

A 

?96 




WRITE (6,143) (I ,YHH ( I ) . 7HH ( 1 ) «PPP ( 1 1 . 1 = 1 ,NTOTT ) 

A 

297 




WRITE (6,144) 

A 

290 




WRITE (6,145) ( I , SUM ( I ) , I = 1 , MOTT ) 

A 

?99 

300 

C 



A 

300 


C 


READ, WRITE CL SREF , ETC 

A 

301 


c 



A 

30? 




REAR (5,136) CL DE S , SREF 

A 

303 




PSAVF. = 6(1) 

A 

304 

305 



IF (N| | INE ,E0.1 ) 00 TO 45 

A 

305 




DO 44 1=?,NLLINF 

A 

306 




HTFM r P ( I ) 

A 

307 



A 4 

6SAVE = AMAX1 (PSAVE ,HTEM) 

A 

308 



45 

COMTIMIlE 

A 

309 

310 



ARAT = R5AVF»PSAVE /SREF 

A 

310 




SHIN = SNM(l) 

A 

311 




00 46 I = ?,MOTT 

A 

312 




STEM = SNN(l) 

A 

313 



46 

SMIN = AMIE'l (SPIN, STEM) 

A 

314 

315 



TOl = 5,F-05»SMN«NSPT (?) 

A 

315 




TOL? = TOL 

A 

316 




TOL 3 = 0.005*SMIN 

A 

317 




WRITE (6,157) TOL? 

A 

310 


c 



A 

319 

320 

c 


f MO OF NEPAL OF OM CALCS 

A 

320 


c 



A 

321 



47 

CON T I MIlE 

A 

322 




IF (lORAG.FO.?) 00 to 5? 

A 

323 




WRITE (6,146) CLOF.S, SREF, RSAVE, ARAT 

A 

324 

3?5 



WRITE (6,145) 

A 

325 




IOPAG = I OP AG* 1 

A 

326 


c 


SET OP ALL A FOR MONK CRITFRION OPTIMIZATION 

A 

327 




00 45 I=1,NT0TT 

A 

328 




00 4 H J=] , N T 0 T T 

A 

329 

330 



CALL WC A| ( 1 , J.MT0TT,YHH,7HH,PPP,A0PT (I ,J) ) 

A 

330 




AOPTd.J) = AOPT ( I ,J) »?./SRtF 

A 

331 



40 

CONT INOF 

A 

332 



40 

COMTIMOE 

A 

333 




IL = NTOTT ♦ 1 

A 

334 

335 



00 50 1 = 1 , MTOT T 

A 

335 




PSIM(I,1> = COS (PRR ( I ) ) 

A 

336 
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34 0 


345 


350 


355 


360 


365 


370 


375 


C 


C 


C 

C 

C 


50 


51 


52 


53 


54 


AOPT (I ,11. J = COS (RRP ( I ) ) 

COM INI IF 

DO 5) 1 = 1, MOTT 

DO 51 .1=1, NTOTT 

AS I M ( J , J ) = AOPT (I, J) 

COMTIHItF 
III = MOTT 

ABOVE 111 IS TEMPORARY SET ER NTOTT 
GO TO 77 
CONTI HUE 

SFT UP ALL A THROUGH F FOR ALL l.J TO DO DIRECT OPTIMIZATION 
I DRAG = I DRAG* 1 
IF lNLl.INE.MF.il GO TO ? 

NTOT2 = NTOTT*? 

NTOT3 = NT OTT * 3 
00 53 I=l,MTOT3 
DO 53 J= 1 , NT OT 3 
T) (l.J) = 0. 

T?( I , J) = 0. 

T3II.JI = 0. 

T4(1,J) = 0. 

T 5 ( I , J ) = 0. 

T 6 ( I , J ) = 0. 

CONTINUE 

00 54 1=1, NTOTT 


II = 1*1 


DO 54 J= 1 .NTOTT 


J1 = J*1 
CALL DRACAL 
T 1 ( 1 1 , J I ) = 
T? ( II ,J1 > = 
T 3 ( 1 1 , J 1 ) = 
T4 ( II ,J1 ) = 

T5(I1,.I1) = 
TF ( II , J1 ) = 
CONTINUE 


( I , J, YMH , 7HM.PPP , A INT I 
AINT(l) 

A I NT ( ? ) 

AINT (3) 

A 1 NT ( 4 ) 

A 1 N T ( 5 ) 

A I N T ( 6 ) 


DIRECT OPTIMIZATION MOOS 

00 56 1=1, NTOTT 

II = 1*1 

DO 56 .1=1 .NTOTT 


A 337 
A 336 
A 339 
A 340 
A 341 
A 342 
A 343 
A 344 
A 345 
A 346 
A 347 
A 348 
A 349 
A 350 
A 351 
A 352 
A 353 
A 354 
A 355 
A 356 
A 357 
A 35B 
A 359 
A 360 
A 361 
A 362 
A 363 
A 364 
A 365 
A 366 
A 367 
A 368 
A 369 
A 370 
A 371 
A 372 
A 373 
A 374 
A 375 
A 376 
A 377 
A 378 


Co 

--4 


PROGRAM ONWASH 



J 1 = J * 1 

A 

379 

3R0 

AOPT (I , J) = T3( 11 ,J1 1-T4 ( 11 ,J1 |-T5( 11 ,J1> ♦T6(TI . J1 ) *T3 < I 1 -1 » JU -T4 

A 

380 


1(11-1, J1>*TS(I1-1.J11-T6(I1-1. Jit. T3M1-1 , Jl-1) +T4 (11-1 • Jl-1) «T5 (I 

A 

381 


21-1, J1-1>*TA(I1-1.J1-1»*T3C11.J 1-11* T* Ml ,J1-1 )-T5(Il,Jl-l)-T6CIl , 

A 

382 


3 J 1 - 1 1 ♦ 1 ,5*SMM( I )« (1 1 (1 1 .J11-T2M 1 ,J1> 1*1 ,5»SNN ( 1 1 « <T1 (Il,Jl-n«T2( 

A 

383 


411 , Jl -1 ) 1 

A 

384 

385 

IF (I.FQ.l) GO TO 55 

A 

385 


AOPT ( T , J) = AOPTM.JI *.5*SNN(I-1 )«M1 M 1-1 .J11-T2M1-) ,J1I*T1(I1-1 

A 

386 


1 ,J1-1 1 *T2 ( I 1-1 , Jl-1) 1 

A 

387 


55 FONT I Milt 

A 

388 


AOPT ( T , J 1 = AOPT 1 1 • J ) *0 ,25 

A 

389 

300 

56 CONTI Ml IF 

A 

390 


PO b4 I = l,f"TOTT 

A 

391 


11 = 1*1 

A 

392 


IF (I.F0.1) GO TO GO 

A 

393 


SCON = 0 ,5* ( SMN Ml* SNN ( I- 1 1 1 

A 

394 

395 

TFMP = 0. 

A 

395 


00 59 J=1 , N T 0 T T 

A 

396 


Jl = J«1 

A 

397 


TEMP = , 5 * SNN ( I -1 ) <* (T 1 ( 1 1 , Jl ) -T2 ( I 1 , Jl ) *T1 ( 1 1 , Jl-1) *T2( 1 1 , Jl-1 1) 

A 

398 


IF (I.EQ.NTOTT) GO TO 58 

A 

399 

400 

IS = 1*1 

A 

400 


no 57 IP=IS,NTOTT 

A 

401 


IP1 = IF * 1 

A 

402 


57 TEMP = TEMP*SCON» (T 1 ( IP1 ,J1 )-T2( IP1 ,J1 1 *T1 ( IP1 , Jl-1 ) *T?( IP1 , Jl-1 1) 

A 

403 


58 CONTINUE 

A 

404 

4 05 

AOPTII.JI = AOPT ( 1 , J) +TEMP 

A 

405 


AOPTM.JI = AOPT ( 1 «J) /SREF 

A 

406 


59 CONTINUE 

A 

407 


GO TO 53 

A 

408 


60 CONTINUE 

A 

409 

4 1 0 

SCON = ,5«SNN(1I 

A 

410 


00 62 J=l,NTOTT 

A 

411 


Jl = J*1 

A 

412 


TEMP = 0. 

A 

413 


00 fcl IP=?,MTOTT 

A 

414 

415 

IP1 = TP ♦ 1 

A 

415 


TEMP = TEMP*SC0N» (T1 ( 1P1 , Jl 1-T2 ( IP1 ,J1) *T1 (IP1 ,J1-1 1 *T2( IP1 ,J1-1 1 1 

A 

416 


61 CONTINUE 

A 

417 


AOPTM.JI = AOPT (I , J) *TFMP 

A 

418 


AOPTM.JI = AOPT (I ,J) /SREF 

A 

419 

420 

6? CONTINUE 

A 

420 



PROGRAM DNWASH 


PAGE 
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63 

CONTINUE 

A 

421 


64 

continue 

A 

422 



00 65 1=) .NTOTT 

A 

423 



00 65 J=l, NTOTT 

A 

424 

425 

65 

T 1 I l , J) = AOPT I J. I ) 

A 

425 



DO 66 1=1, NTOTT 

A 

426 



00 66 J=l, NTOTT 

A 

427 



AOPTII.J) = ?.« (AOPT ( I . J) *T 1 ( I • J)) 

A 

428 


66 

CONTINUE 

A 

429 

430 


IL = NTOTT. I 

A 

430 



III = II. 

A 

431 



IL2 = IL* 1 

A 

432 



*3 

O 

-<4 

II 

A 

433 



AOPTII.IL?) = 0. 

A 

434 

4 35 

67 

CONTINUE 

A 

435 



A0PTIIL1.il?) = CIOFS 

A 

436 



AOPT { IL I ♦ 1 ) = P*COS (PPP (1)1 »SNN ! 1 )•«?/( 3 . »SRE F ) 

A 

437 



DO 61! 1 = ?, NTOTT 

A 

438 



AOPT I IL 1.1) - (2<*C05 (PPP ( I ) )«SNN(I)«»2*C0S(PPP(I-1) ) »SNN ( I - 1 ) «•?) * 

A 

439 

440 


1 4 . / 1 3 , <*SRE F ) 

A 

440 


6fl 

CONTINUE 

A 

441 



SOM* = 0. 

A 

442 



00 69 1=?, NTOTT 

A 

443 



SOM* = SHE’* *COS. (PPP (II) *SNN ( I ) 

A 

444 

445 

69 

CONTINUE 

A 

445 



AOPT ( IL 1 , 1 ) = AOPIIII 1,1)*|4./SREF)«SNN(1)»SUMX 

A 

446 



00 7? 1=?, NTOTT 

A 

447 



SON* = COS(PPPIl) )«SNN(I)«SNN(I-1I 

A 

448 



11 = 1*1 

A 

449 

450 


IF (I. FO. NTOTT) GO TO 71 

A 

450 



1)0 70 J=I 1 .NTOTT 

A 

451 



SOMX = SOM X * COS (PPP 1 J) 1 *5NN ( J) * (SNN 1 1 ) *SNN ( I -1 ) ) 

A 

452 


70 

CONTINUE 

A 

453 


7] 

CONTINUE 

A 

454 

455 


AOPT(ILl.l) = AOPTIIl l,I)*(4./SRtF)*SlmX 

A 

455 


72 

CONTINUE 

A 

456 


C 


A 

457 



DO 73 1=1, NTOTT 

A 

458 



AOPT(lLl.I) = ? . » AOPT I IL 1 . I ) 

A 

459 

4 GO 

73 

CONTINUE 

A 

460 



00 74 1=1, NTOTT 

A 

461 



A0PTU.IL1) = AOPT ( II. 1 , I ) 

A 

462 


w 

VO 



4 > 

O 


PROP, HAM DMWASH 


PARE 


1? 


74 CONTINUE A 463 

00 75 1 = 1 ,11.1 A 464 

465 CLP ( I ) = AOPT ( ! 1. 1 • I ) A 465 

75 CONTINUE A 466 

WPITF (6,147) (CLP ( T ) . 1=1 , IL 1) A 467 

AOPT ( II. 1 r II 11=0. A 468 

DO 76 1=1, 1LI A 469 

470 00 76 J= 1 , 1 L 1 A 470 

AS1M(I,J» = AOPT(I.J) A 471 

t'Sl N |I,1) = AOPT ( I , IL?) A 47? 

76 CONTINUE A 473 

WRITE ( 6 , 146 ) A 474 

475 77 CONTINUE A *75 

CALL SIMEO (ASIM, IL 1 , HS IM , 1 ,DE T , I P I VOT , 51 , I SCALE. ) A 476 

IT = IL l ♦ l A 477 

00 70 1=1, IL1 A 478 

78 AOPT(I.IT) = PS!M(1,1) A 479 

480 C SET I L 1 PACK TO NTOT T ♦ 1 A 480 

IL1 = N T 0 T T ♦ 1 A 401 

IF (InPAC-.ME.31 GO TO 81 A 48? 

Cl.CHK = 0. A 483 

no 79 I=l,NTOTT A 484 

485 AOPT ( I , IL?) = BSIM(l,l) A 485 

CLCHK = CLCHK*CLP ( I ) «A0PT ( I , IL?) A 486 

79 CONTINtlF A 487 

VIRITE (6,150) CLCHK A 488 

DO fill 1 = 1,11 A 409 

490 80 WRITE (6,15?) 1 , AOPT [ I , I L ? ) A 490 

81 CONTINUE A 491 

WRITE (6.151) A 49? 

C A 493 

C CALCULATE HOUND CIRCULATIONS AND A 494 

495 C CALCULATE CMTPl N IIOWMWASH DIVIDED R V FRFESTREAM U A 495 

C A 496 

NTOT 1 = NTOTT ♦ I A 49T 

SUMGAM =0. A 49B 

AOPT (NTOT1 , IL) =0. A 499 

500 IF (IPRAG.NE.3) GO TO 83 A 500 

00 82 1=1, IL A 501 

fl? AOPT ( 1 , IL 1 = Appling?) A 50? 

AOPMIL.IL) =0. A 503 

63 CONTINUE A 504 
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505 


1)0 84 I - 1 , NTOT 1 

A 

505 



CALL GAMCAL ( 1 , 1L ,SNN , AOP 1 , BGAMO 1 

A 

506 



BGAMII] = BGAMO 

A 

507 


64 

CONT I Ml IF 

A 

508 



IE (NIL IMF. ,F0.1» GO TO 90 

A 

509 

510 


ML 1 = NLL INF-1 

A 

510 



00 85 J=l, Nil 

A 

511 



11 = Ml t INF - I • 1 

A 

512 



11 = ! F ♦ 1 

A 

513 



JF = NSPT 011-1 

A 

514 

515 


GSllM(I) = HGAM ( JF) »AOPT ( JF , 1L I »SNN ( JE) 

A 

515 


85 

CONT I Milt 

A 

516 



DO 89 1 = 1, ML 1 

A 

517 



11 = Ml LFMF-I*? 

A 

518 



JF = NSPT ( F I ) - I 

A 

519 

520 


YT = YHH ( Jf | »5MN ( JF ) »COS (PPP ( JE )) 

A 

520 



IF (A8S1Y1 1 .1 T .O.OUOl) GO TO 89 

A 

521 



JF = N5PT l 2) - 1 

A 

522 



00 86 J=1 , JF 

A 

523 



JS = J 

A 

524 

525 


IF (YT.LT.YHH(J)) GO TO 87 

A 

525 


86 

CONTINUE 

A 

526 


87 

CONTINUE 

A 

527 



DO 88 J=JS,Jf 

A 

52B 



RGAM(J) = RGAN ( J I » GSUM ( I ) 

A 

529 

530 

ea 

CONTINUE 

A 

530 


85 

CONTINUE 

A 

531 


90 

CONT INI IF 

A 

532 



SIIMGAM = 0. 

A 

533 



DO 94 I=l,NTOTT 

A 

534 

535 


DO 91 ,)=1 ,NLl. INF 

A 

535 



JJ = J*1 

A 

536 



JChK = MSPT(JJ)-1 

A 

537 



IF (I ,FO. JCMK ) GO TO 9? 

A 

538 


91 

CONTINUE 

A 

539 

540 


SIIMGAM = SIIMGAM ,COS( PPP ( I) ) * ( SNN O ) »«2» I AOPT { I « 1 * ID ♦2,*A0PT ( I , IL) 

A 

540 



I 1 /3, « SNN ( I ) ef'GAM (III 

A 

541 



GO TO 93 

A 

542 


92 

CONTlUUf- 

A 

543 



S1IMGAM = SIIMGAM, COS (PPP (1 1) » (SNN( 1 ) *«2"2.«AOPT 1 1 , !L> /3.*SNN 1 1 ) *BGA 

A 

544 

545 


1H ( I > » 

A 

545 


93 

CONTINUE 

A 

546 



PROGRAM OMWASH 


94 CONTINUE A 647 

IF (NLI INE.MF.11 GO TO 95 A 548 

JE = NSPT (?) -1 A 549 

550 95 CONTINUE A 550 

RGAM(NTOTl) = ROAM ( JE ) ♦ AOP T ( JE * IL) *SNN ( JE ) A 551 

VlPHU = CLOE5»SRFF/(H.«SUMGAM) A 552 

C A 553 

C PtSCALF OPT SHFI) SMFET STRENGTHS A 554 

555 C TO 0E DIVIDED RY I) INSTEAD OF W A 555 

C A 556 

IF (IDRA0.F0.3) fid TO 97 A 557 

DO 96 1=1 ,MT0T1 A 556 

BfiAM ( I ) = HOAR ( I ) “WPfUl A 559 

560 96 AOPTCUIL) = AOP T ( 1 , 1 L I «WOHll A 560 

97 CONTINUE A 561 

C A 562 

C HfSCAlF HOUND ♦ SHED STRENGTHS TO CALCUl ATE AVG ♦ NON-DIM VALUES A 563 

C A 564 

565 C AOPT(I.IL)=OPT SHED SHEET STRENGTHS A 565 

C HGAM ( I ) =OPT POUND CIPC VALUES A 566 

C CDPAG ( T 1 =OPT NON DIM HOUND CIRCS A 567 

C GAM ( I 1 = opt NON DIM SHEO SHEET VALUES A 568 

C A 569 

570 SUBGAM =0. A 570 

SUHGAM =0. A 571 

MS = NSPT (2)-] A 572 

DO 100 I-l.NS A 573 

TFMP = AOP T (1*1,11 » A 574 

575 DO 99 J= 1 ,NLL INF A 575 

JJ = J*1 A 576 

jrHK = NSPT I JJ) -1 A 577 

IF (I.f.F.JCHK) GO TO 9fl A 578 

TEMP = 0 . A 579 

580 9p CONTINUE A 580 

99 CONTINUE A 5B1 

SUKGAM = SUMGAh*SMN ( I ) « (TEMP* AOPT ( I , IL) ) A 582 

SUHGAm = SUGGAM«2. »SNN( I ) # (HGAM ( 1 ) ♦ (SNN ( I ) /12. ) • ( 4 , <• TFMP*B . *A0PT ( I A 583 

I , IL ) ) ) A 584 

585 ion CONTINUE A 585 

SUMGAM = SUMf-AM/OTUU A 586 

SU 6 GAM = SUHfiAM/OTOn) A 587 

NS = NSPT (?)-? A 588 
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M51 = N'5+l 

A 

5B9 

590 


1)0 103 1 = 1. PS 

A 

590 



IF (ICNTPt .FO.M GO TO 101 

A 

591 



CDPAO (1*1) = PG AM ( 1 ♦ 1 ) /SUPGAP 

A 

592 



GO TO 10? 

A 

593 


101 

CONTINUE 

A 

59* 

595 


COP AG ( I + 1 1 = PGAM (1*1) /HGAM (NS1 1 

A 

595 


10? 

CONT IMIlF 

A 

596 


103 

CONTINUE 

A 

597 



Jl = N5P T ( ? ) - 1 

A 

598 



COPAG(NTOTl) = HGAP < JL) tAOPT ( JL. IL) »SNN ( JL) 

A 

599 

600 


COPAG(NTOTl) = C0PAG(NT0T1)/SUH6AH 

A 

600 



00 104 1= l .NTOTT 

A 

601 



GAP (I) = AOPT(I.IL) /SUPGAP 

A 

602 


1 04 

COMTINOt 

A 

603 



COPAG(l) = 0. 

A 

60* 

605 


GAP (NTOT 1 ) = 0. 

A 

605 



IF (NI.L1NE .f 0.1 ) GO TO 109 

A 

606 



00 10P I =2 .Nl LINE 

A 

607 



J5 = N5PT ( I ) 

A 

608 



II = 1*1 

A 

609 

610 


JF = NSPT ( T I ) -1 

A 

610 



500 GAP = 0. 

A 

611 



00 106 J=J5,JF 

A 

612 



TFPP = AOPT (J.l . II ) 

A 

613 



IF (J.NE.Jt) GO TO 105 

A 

61* 

615 


TFPP = 0. 

A 

615 


105 

CONTINUE 

A 

616 



SIIHGAP = SUPGAP *?.»SNN ( J) • (PGAM ( J) ♦ (SNN(J) /!?.)• ( * . *TEMP* 8 . • AOPT ( J 

A 

617 


l.ILH) 

A 

618 


106 

CONTINUE 

A 

619 

6?0 


SUPGAP = SIIPf.AP/OTO ( I ) 

A 

620 



00 107 J= JS • JF 

A 

621 



CnRAG(J) = PGAP (J)/SI)BGAM 

A 

622 


107 

CONT1NIIF 

A 

623 


10P 

CONTINUE 

A 

62* 

6?5 

109 

CONTINUE 

A 

625 



no 11? 1=1. NTOT] 

A 

626 



IF (I.F0.NT0T1) GO TO 1 10 

A 

627 



FTA = ?.» (-YE'P ( I ) *SNN(I ) *COS(PPP( I )) ) /HSAvE 

A 

6?8 



GO TO 111 

A 

629 

630 

110 

CONTINUE 

A 

630 


■P- 

w 



■> 
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635 


64 0 


645 


650 


655 


660 


665 


670 


C 

C 

C 


FTA = 0. 

Ill CONTINUE 

WRITE (6.15?) I ,PGAM ( I ) .CDRAG ( 1 ) . AOPT ( I . IL) .GAM ( I ) , FTA 
] 1? CONTINIIF 

INOUCFO DRAG CALC USING ANALYTICAL INT OT WN TIMES ROUND CIRC 


CD = 0. 

NI L = NLL 1 NT ♦ 1 

DO 11R I=l,MTOTT 

DO 113 K=?,MLl 

ICHK r NSPT ( N J - 1 

If (ICHK. TO. I) GO TO 114 

113 CONTI Milt 

XI = RGAM ( I ) ♦ ,?5»SNN ( I ) » ( AOPT ( I ♦ I • IL) O.^AOPT ( I » IL) » 

Y I = .5® (AOPT ( I ♦ 1 « IL ) ♦ AOPT ( I • IL ) ) 

71 = .5® (AOPT (1*1 .ID-AOPT ( I .ID ) 

GO TO 115 

114 CONTINUE 

XI = RGAM ( I ) ♦,75«SUN ( I ) «AOPT ( I . IL ) 

Yl = . 5® AOPT ( I . IL ) 

Z! = -Y I 

115 CONTI NUT 

DO 119 J=1 .NTOTT 
DO 116 K = ?,NI l 
JCHK = NSPT (M ) - 1 
IT (J.FQ.JCHK) GO TO 117 

116 CONTINIIF 

YJ = .5® (AOPT ( J* 1 . IL) *A()PT ( J. IL ) ) 

7J = .5» (AOPT (J. 1 , II l-AOPT (J. IL) ) 

GO TO IIP 

117 CONTINIIF 

YJ = 0 .5® AOPT (J.II ) 

7J = -YJ 

118 CONTINIIF 

CALL DPACAL I I . J , YHH , 7HH ,PPP . A INT ) 

CD 1 = XI *YJ»«INT ( 1 ) ♦» I *7 J« A INT (2) ♦YI»YJ*AINT (3) «YI«7J®AINT (4) »ZI«Y 


1 J*AINT(5) ♦7I*7J»AINT(6) 
CD = CD *C0 1 
HR CONTINUE 

CD = CD® 2 • 7 5REF 
IF ( IDRAG.E0.3) 


A 631 
A 63? 
A 633 
A 634 
A 635 
A 636 
A 637 
A 638 
A 639 
A 640 
A 641 
A 64? 
A 643 
A 644 
A 645 
A 646 
A 647 
A 648 
A 649 
A 650 
A 651 
A 65? 
A 653 
A 654 
A 655 
A 656 
A 657 
A 658 
A 659 
A 660 
A 661 
A 66? 
A 663 
A 664 
A 665 
A 666 
A 667 
A 668 
A 669 
A 670 
A 671 
A 67? 




GO TO 1 20 
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WRITE (6,154) CO 

A 

673 



GO TO 121 

A 

674 

675 

1 20 

CONT INIlE 

A 

675 



VIPITI (6,153) CO 

A 

676 


121 

CONT 1 Hilt 

A 

677 



DIDfAl = SRFFofl OfSe«2/(Pl«(HSAVt)in»2) 

A 

678 



IF (IOPAG.FQ.3) GO TO 122 

A 

679 

680 


WRITE (6,155) Pint At , WOHO 

A 

680 


122 

CONT I Ml It 

A 

601 



OEFF = 1)11)1 AL/CU 

A 

682 



WRITE (6.156) PFEE 

A 

603 



ISTOP = ISlOP.l 

A 

684 

685 


PO 123 I=1,NT0T1 

A 

685 


123 

GAM (I) = AOPT ( I , II ) 

A 

606 



C-0 TO 3 

A 

687 


124 

CONI I HUE 

A 

688 



CDAPP = 0. 

A 

689 

690 


PO 120 I=l,NTOTT 

A 

690 



DO 125 J=1 rW.I IMF 

A 

691 



JJ = J,1 

A 

692 



JChK s NSFT(JJ)-) 

A 

693 



IF (l.FO.JCHK) GO TO 126 

A 

694 

695 

125 

CONTINUE 

A 

695 



CDAPP = C|IAPP«2 . «SNM ( 1 1 «»NWSH ( I ) » (8GAM ( I ) ♦ , 25« ( AOPT (1*1 , ID ♦3,»A0P 

A 

696 


IT (I ,11 ) )«SNN<1 1 1 

A 

697 



GO TO 127 

A 

698 


126 

CONTINUE 

A 

699 

700 


CDAPP = CDAPP + 2 , B 6NN ( 1 ) <M)NWSH ( I)« (HGAM ( I ) *.75»A0PT ( I ,IL1*SNN(I) ) 

A 

700 


127 

CONT IMi/E 

A 

701 


128 

CONTINUE 

A 

702 



CDAPP = 4 .«CPAPP/5REF 

A 

703 



WRITE (6,158) CDAPP 

A 

704 

705 


IF (10RAG.E0.3) GO TO 2 

A 

705 



•a- 

C 

C 

A 

706 


1 29 

CONTINUE 

A 

707 




A 

708 


130 

FORMAT ( 1H1 ) 

A 

709 

710 

131 

F OP M A T ( 5 1 5 ) 

A 

710 


132 

FORMAT (25X.3F15.5) 

A 

711 


133 

format (//33X.1 HI, 7X,8HnOv;NWASH,7A,10HW/CO5 (PHI) /) 

A 

712 


134 

FORMAT ( 3 OX . I5.2F 15.5) 

A 

713 


135 

FOPMAT (/30X,?2HGFNFRAL INPUT GEOHETRT//) 

A 

714 


-P> 

t_n 



-O' 

a) 
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T 1 5 Uf. FOPMAT fbFIO.O) A 715 

1 37 FORMAT ( 1015) A 716 

1 3P FOPMAT(30X..37HT()TAL PLANFOPM PEP1PHEPAI. LENGTH= , F 1 5 .5/ ) A 717 

139 FORMAT (/30X , 1 7HCIRCULAR ARC WING) A 718 

1 A 0 FORMAT (/30X.21HFOUAL SFGMFnT SPACING) A 719 

720 1 A 1 FORMAT (/30X.22HCOSIMF St GHENT SPACING) A 720 

14? FOPMAT </2AX .PHSFRMT NO. 8X , 1 HY • 15X . 1 HZ » 1 3X » 3HPHI / ) A 721 

1A3 FORMAT (25X. 1S.3F15.6) A 722 

1 A A FORMAT ( /?7X , 1HI ,HX .6HSNN ( I )/ ) A 723 

1 A5 FOPMAT (25X , 15, F 12.5) A 72A 

725 1A6 FOPMAT!/// ?5X , 20HOFS I GN LIFT COEF. = , F 1 0 . 5// , 25x , 22HW ING REFEREN A 725 

ICE AREA = *F 1 0.5//* ?5x , 16HREF WING SPAN = , F10.5//.25X, A 726 

215HASPFCT RATIO = ,F10.5/) A 727 

1 A 7 FORMAT (/IX ,) OF 1 0.3) A 728 

1 A8 FORMAT ( ///?5x , s AHO 1 Rf CT OPTIMIZATION USING ANALYTICAL EXPRESSION F A 729 
730 1 OR Cl)//) A 730 

1A9 FORMAT (///2SX,37HOPTIMUM LOADING USING HUNKS CRITERION//) A 731 

150 FORMAT (//23X.AFHLIFT COFFF CALCULATED FROM CLP AND SOLVED GAM=, A 732 

IE 1 3 . 5/ ) A 733 

151 FORMAT <//??X.5PSFGMt,3X,10HGOUND CIRC,2X,BHDGAM/AVE,AX»10HSHED STR A 73A 

735 1H.SX.7UGAM/AVE , Rx , 3|iF T A// ) A 735 

152 FORMAT <?0X,I5, C F13.R) A 736 

153 FOPMAT !///?SX ,( 7HCO CALCULATE 0 USING SUB DRACAL AND LOADS FROM DIR A 737 

1ECT OPTIMIZATION = .F1S.S//) A 738 

ISA FORMAT ( ///25X ,<4HCD CALCULATED USING SUB DRACAL AND OPTIM LOADS US A 739 
740 11HG HUNK CRIT =,F1S.5//1 A 740 

155 FOPMAT (25X, 19HCO FOR FLAT WING = ,E 1 5 . 5/25X . 35HRAT 1 0 OF ZERO DINED A 741 

1PAI nnwMWSH/ll = .FIS. 5/) A 742 

156 FOPMAT (//25X,49MIM)UCED DRAG EFFICIENCY FOR WINGS OF EQUAL SPAN = A 743 

1 .F IS. 6/) A 744 

157 format (/25X.14HTOI IN SNT AN =,E15.5/) A 745 

158 F0RMAT</25X,?RPAPPROx CD USING SOLVED ROUND/25X ,4?HC IRCULAT IONS AN A 746 

ID W A SHFS AT SFP MIDPOINTS =,E13.5/) A 747 

END A 748 


745 


noon 
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I 


in 


15 


20 


25 


SUBHOIIT INF GAMCAL ( 1 < 11. » SNN i AOPT » HOAMO ) 

CAl.CMI ATE BOUND C 1 POOL A T I ON AT Lt F T EM) PT OF SEGMENT l.BGAM 
AOPT LAST cm HUN contains APB ay OF OPTIMIZED SHED SHEET STRENGTHS 

COMMON /FEN/ NSP T ( 1 0 ) t Nl | |NE 
DIMENSION SNN(l), AOP T (52* 1 ) 

DC-AMO = 0. 

no i k = i.ni i ine 
KK = N| LINE -Ml 
KCHK = NSPT (KK) 

IF (I .FO.KCUK) GO TO A 
IF (I. OF. KCHK) GO TO 2 

1 CONTINUE 

2 CONTINUE 
KST = KChK 

nOAMfl = AOUT (KST , IL ) «SNN (KST) ♦AOPT (1 . IL) »SNN( 1-1 ) 

KCHK1 = K CHK ♦ 1 
IF (I.FO.KCHK1) GO TO 5 
TM = 1-1 

IP = KCHK ♦ 1 

DO 3 J=1P,IH 

6GAM0 = GGAMO^AOPT ( J » IL)* (SNNIJ-l ) ♦ SNN f J) ) 

3 COUTIMUF 
GO TO 5 

A 8GAM0 = n, 

5 CONTINUE 
PFTIIBN 
FNO 


H 1 
B 2 
B 3 
B 4 
B 5 
B 6 
B 7 
B B 
B 9 
B 10 
B 11 
B 12 
B 13 
B 14 
B 15 
B 16 
B 17 
B 18 
B 19 
B 20 
B 21 
B 22 
B 23 
B 24 
B 25 
B 26 
B 27 
B 28 
B 29 


■C- 



•p, 

00 


SUBROUTINE WCAt 


PARE 


1 


] 


5 


10 




an 


25 


30 


35 


AO 


C 

c 

c 

c 

c 


SUBROUTINE WCAL ( I • J.NTOT , yHH f 2HH ,PPP « AA A A 1 C 

C 

CAI.CIHATE MATRIX COEFFICIENT FOR DRAG OPTIMIZATION USING MONK CPIT C 
FINDS COfF OF I FH DOWNWASH DUE TO JTH SHED SHEET STRENGTH C 

IF, FINDS COfF MULTIPLYING STRtNGIH J IN EQUATION I C 

C 

DIMENSION YHH(l), 7HH(1I, PPP(l) C 

CORDON /SEC,/ SNN (51) C 

CONDOM /FEN/ NSPT ( 10) «NLI. INF C 

INTFGFR R C 

PI = 4.»ATAh(l.) C 

A AAA = 0. C 

ICNTRL =0 C 

P - I C 

K = J C 

1 CONTINUE C 

CALI CC At (P,K , YHH, ?HH,PRP,SNN (K ) , AA ,RB,DD,FF ,GG,EE « AJ, AK,RR,TT ,UU C 

I.WW) C 

CALL CONCAL ( A A , HH , FF , GG , SNN (K ) , A , B ,C ,0 ,E , G • H J, BK , RL , RM * BN , BO ,BP , 1 C 
1) C 

IF (RR.FO.O.) GO TO 2 C 

EXPR = ? .» ( ATAN? (C, AHSIRRI) -ATAN2 <D,AHS (PRD) /(ABS(PR) 1 C 

GO TO 3 C 

2 CONTINUE C 

EXPR = 2./D-P./C C 

3 CONTINUE C 

FII.OG = Al.OG (F/G ) C 

A 1 PK = (A»FXPR*.S»HH«RLOG)/(?.«PI) C 

A 3PK = (HL»Expp«2.«HP»SNNtK) ♦B0«RL0G)/(2.»PI*SNN(K) » C 

CA|L CONCAl |I)D,FF , A J, AK , SNN (K ) . A ,H ,C ,D.F ,G.B J, BX ,BL .BM»BN, BOt BP , 2 C 

1 ) C 

IF (Ull.bO.O.) GO TO A C 

F.XPP = 2. » ( AT AN2 (C. APS (III))) -ATAN2 (D,ADS IUU) ) ) / I ABSIUll) ) C 

GO TO S C 

4 CONTINUE C 

EXPR = d./V-P./C C 

5 CONTINUE C 

PL 00 = ALOG (F/G I C 

A?PK = - ( A *FXPR ♦ , 5*E F *RLOG) / ( 2, *P 1 1 C 

A4PK = -IK «EXPP*2."EE*5NN(K) «B0«RL0G1/ t?.»PI*SNN|Kl > C 

IF (ICNTRL. FO.?) GO TO 7 C 

A A AA = ( A IPX « A?P|< ) » .SO- ( A3PK, AmPK ) » .50 C 


1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

n 

12 

13 

14 

15 

16 
IT 
IB 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 


SlIRROUT INF. WCAL 


PAGE 


45 


50 


ICNTRI. = ? 

00 6 l. = l .NILTNF 
JCHK = NSPTd.) 

IF (K.FO.JCHK) GO 10 H 

6 CONTINUE 
K = J-] 

GO TO I 

7 CONTINUE 

AAAA = AAAA+O.S* (A|RK+A?PK*A3HK+A4PKl 
fl CONTINtlf 
RETURN 
END 


C 43 
C 44 
C 45 
C 46 
C 47 
C 4B 
C 49 
C 50 
C 51 
C 52 
C 53 
C 54 
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SUBROUTINE CCAL 


PAGE 1 


] 


10 


15 


20 


25 


30 


C 

C 

C 

C 

C 

C 


SUBROUTINE CCAl ( I ,.J, YHH, ZHH.PPP.Sf AA,flG,ODiFF,GG,EE, JJ,KK,RR,TT,U 
1U.HW) 

SURROIITINF CCAL 

calculates geometrical constants needed in evaluation of integrals 

FOR VARYING I AND J VALDES 
REAL JJ.fcK 

DIMENSION YMh(l). 7HH(1), PPP(l) 

DY1J = YHlUn-YHHIJ) 

07 1 J = ZHH»I)-7HH(J) 

COI = COS(PPPII) ) 

SI I = SIN (PPP ( I ) ) 

COJ = COSIPPPIJI) 

SIJ = SIN (PPP (J! ) 

AA = DY1J»C0I*P71J*SI I 

rr = -costpppui-ppptm 

Ff = (OYI J«C0J*I)7I J°SI J) 

GG = OYl J»PYl J*D7l J»DZ1J 

DYIJP = YHM ( T ) ♦ YMH ( J J 

OZIJP = 0Z1J 

DD = DYlJPoCOI ♦OZlJPoSII 

Ff = COSIPPP(J) + PPP (1)1 

JJ = 2.<MOYIJP»COJ-DZIJP»SIJ) 

KK = DyI JP<tnYIJP*D7IJP«DZIJP 
PR = (DY I J»S I J-DZT J*COJ) 

TT = 2»SIN (PPP ( J) -PPP ( I )) 

11(1 = 2.»<DYIJP«SI0*0ZIJP«C0J) 

HW = 2. <*SIN(PPP ( J) *PPP < I ) ) 

RETURN 

END 


D 1 
0 2 
0 3 

D A 
D 5 
D 6 
D T 
0 0 
D 9 
D 10 
D 11 
D 12 
D 13 
D 14 
D 15 
D 16 
0 17 

0 10 
D 19 
D 20 
D 21 
D 22 
0 23 

D 24 
D 25 
D 26 
0 27 

D 20 
D 29 
D 30 
D 31 
D 32 
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1 SUBROUTINE CONCAL ( A A , 86 . FF , GG , S . A ,B ,C ,0 , F , 0 , J, K ,L , M , N, 0,P, ICNTRL ) E 1 

E 2 

SUBROUTINE CONCAL E 3 

E A 

CALCULATES GFOPF TR I CAL CONSTANTS NEEDED IN EVALUATION OF INTEGRALS E 5 

E 

RFAL J»K #L t M »N E 7 

A = AA-0.5*BP«FF E B 

B = l.-HH*HR E 9 

10 C = FF*2.*S E 10 

D = FF-2 . *S E H 

F = S*S«S«FF*GG E 12 

G = S»S-S«FF«GC- E 13 

J = 2.«<AA*S«PP) E 14 

15 K = ?.«(AA-S»HH) E 15 

L = 0.5* <fiB»FF*FF-AA»FF-2.«Hb«GG) E 16 

M = 0.5* <-FF-6.*AA»RB*A.»FF*BB*BB) E 17 

N = 2.*<HB»HH-1 .)«RB E 18 

0 = 0.5* ( AA-FF »HR ) E 19 

20 P = 0.5* ( 1 »-?.*HB*BH) E 20 

RETURN E 21 

END E 22 
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SUBROUTINE SNTAN 
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I 


5 


in 


15 


20 


25 


30 


35 


AO 


C 

C 

C 
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SUBROUTINE SMTAN < S , C ,BB . PR . T 7 ,RT AN . RST AN , RS2T AN, RS3T AN ,RSAT AN ) 
SUBROUTINE SNTAN 

EVALUATES INTEGRALS OF THE FORM S««N»ATAN ( (C,2»BB*S) / (PR*S«TT) ) 
ALL OIVIPEU BY <RR»S»TT) 

WITH RFSPECT TO S BETWEEN LIMITS OF -S AND S FOR N=0, 1,2,3. 

ATAN PART OF INTEGRANO APPROXIMATED AS A QUADRATIC IN S WHICH IS 
FORCED THROUGH ATAN VAUJFS AT -S.O.ANU S. 

A , C ARF CALCULATED IN SUBROUTINE CONCAL 
HB.RR.TT ARF CALCULATED IN SUBROUTINE CCAL 

RESULTS APE RT AN.RSTAN , RS2 T AN ,RS3T AN 

APPROXIMATE INTEGRAL EVALUATED USING MACS YMA PROGRAM OF MIT PROJ. 
EVALUATION OF INTEGRALS FOR TT=0. REGIND AT AT LABEL 10 
SINGULAR INTEGRALS EVALUATED AT APPROXIMATE ENDPOINTS , 4-SAW AY 
MIDRANGE SINGULARITIES E XCLUDED , ATAN PART OF INT APPROX- 
IMATED AS 2 QUADRATICS 

DIMENSION AA] 13,31 , AA<3), IPIVOT(3) 

COMMON /TELL/ TOL.TOL2 
RPR = RP 

sss = s 

CCC = C 
RT AN = 0. 

RSTAN = 0. 

RS2TAN = 0. 

RS3TAN = 0. 

RSATAN = 0.0 

IF (ABS(TT) .LT,lE-05.AND.ABS(RRX.LT.lE-05) GO TO 7 
IF (TT.EQ.O.O.AND.RP.EQ.O.OI 60 TO 7 
IF (TT.EQ.O.O) GO TO A 

FIRST, CHECK FOR MIDRANGE S I NGUL AR I T I ES ,E XCLUOING ANY FOUND 
SZERO = -RR/TT 

IF (ABS(ABS(SZFRO)-SI ,LT . IF -03 . AND . ABS ( SZERO ) .LE.S) GO TO 3 
IF (S7FRO.LT.O.O.ANO.S7ERO.GT.-S) GO TO 16 
IF (S7ERO.GE. 0.0. AND. SZERO, LT.SJ GO TO 16 
1 CONTINUE 
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F B 
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F 13 
F 1A 
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F 22 
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F 2A 
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F 26 
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F 31 
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F 35 
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F 37 
F 38 
F 39 
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F A1 
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C = ATAN? (CCC,ARS ( RR I ) F *3 

Cl = ,5«<ATAN2< <CCC®2.«BH*S) ♦ ABS (RR®TT*S >) ®ATAN2 MCCC-?.*BB*SI »ARS F AA 

AS 1(RR-TT*S) ) ) -r F AS 

Cl = C1/<S®S) F A6 

C? = (ATAM?( (CCC*?.»PB«S) ,ARS (RR»TT»S)) )/S-C/S-Cl*S F AT 

C F AB 

C INTEGRANO NOW IS <C1«S»S»C?»S*C)»S®*N/(RP*TT*S) ) F A9 

50 C F 50 

CLOGR = ALOG(ABS((PR + TT*S|/(RR-TT*SH) F 51 

CONO = (C®TT®TT-C2*PR*TT+C1®RR*RR)/(TT*®3) F 52 

CONI = <C2*TT-C1®RR)/(TT®TT) F 53 

RTAN * CONO*CLOGR*?.»S»CON1 F 5A 

55 CON2 = <C»TT®TT-C?*RR®TT®C1«RR»RH)/(TT«®3) F 55 

CON3 = <C«RR*TT*TT-C?»RR»RR»TT*C1*RR»®3)/(TT®»A> F 56 

RSTAN = 2 . ®S*CON2-CON3*CLOGR *2 » ®Cl®TT*TT®S*®3/<3*TT # *3) F 57 

CONA = (C®RP«RR*TT®TT-C2*RR*®3»TT+C1*RR»®A)/(TT # ®5) F 5B 

CON5 = A.»C2*TT»®3-A.®C1*RR®TT*®2 F 59 

60 C0N5 = C0N5/(1?®TT*®AI F 60 

C0N6 = (-C»PR»TT»TT*C2*RR»RR»TT-C1«RR»*3) /(TT»*A) F 61 

RS2TAN = CnNA»CL0PR*2.»C0N5»5»«3*2.»5«C0N6 F 62 

CONB = 20.®(C®TT»«A-C2»RR®TT*®3>C1»RR*®2*TT' , *?)/(60»TT®«5) F 63 

C0N9 = (C*RR»«2*TT«»2-C2»RR«»3»TT*C1»RR«»A)/(TT»*5) F 6A 

65 CON 1 0 = C*PR»®3«TT«»?-C2*RR®®A®TT®C1®RR®»5 F 65 

CONI 0 = CONI 0/ < TT»»6) F 66 

RS3TAN = 2 . »CONB*5**3®2 • *C0N9*5-C0N1 0*CLOGR+2A.*C1®TT®*A®S®*5/(60* F 67 

1TT°“5) F 68 

COMA = <C®RR»«A®TT®*2-CP*RR®®5»TT®C1»RR»®6)/'(TT®®7> F 69 

70 CONH = <C2«TT-C1®RP)/(5®TT®*2) F 70 

CONC = (-C»PR®TT»»2«C2*RP«»2*TT-C1»RR«»3)/(3»TT**A) F 71 

CONO = CONC*RP»»2/ (TT«*2/3) F 72 

PSATAN = CONA "CLOGR ♦CONH® 2* 5* *5 ♦CONC«2®S®»3*COND®2*S F 73 

IF (RR.GT.O.) GO TO 2 F 7A 

75 RTAN = -RTAN F 75 

RSTAN = -PSTAN F 76 

PS2TAN = -RS2TAN F 77 

RS3TAN = -RS3TAN F 78 

PSATAN = -PSATAN F 79 

80 2 CONTINUE F 80 

GO TO 13 F 81 

3 CONTINUE F 82 

SAWAY = S-TOL F 83 

S = SAWAY F 8A 


Ui 

w 


Ln 

-(n 


SUBROUTINE SNTAN 


PAGE 


3 


85 


90 


95 


100 


105 


110 


115 


120 
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c 
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GO TO 1 

FOR CASE OF PR NOT ZERO.TT=0.0 

A CONTINUE 
PR = AHS(RP) 

ALNNUM = t?*RH)*®2»S““2*4*C»RB«S*RR“*2*C«*2 
ALNDF.N = (2*RP) «»2»S*“2-4“C*BB«S*PR“*Z*C*»2 
IF (ALNNUM. FO. 0.0. OR. ALNOFN. EQ. 0.01 GO TO 5 
GO TO ft 

5 CONTINUE 
S = S-TOL 

ALNNUM = <2«PP) •«2«S*“2+4»C»BB*S*PR*“2*C**2 
ALNDEN = <2*RR)““?“S**2”4»C*6B»S*RP“*2“C*“2 

6 RATLH = ALOG ( ALNNUM /ALNOEN ) 

TNOIF = ATAN2I (C*2.“RH»S) ,RR)-ATAN2< <C-2.«BB“SI .RR) 

TMSUM = ATAN2( (C*?.“RB“S) .PR) *ATAN2( (C-2,*BB*S) iRR) 

RT AN = -(.25“PP/RB(«PATLN*0 . 5*C*TN0IF/BH 
PTAN = RTAN/RP*S*TNSUM/RR 

RSTAN = 0.5“ (S“S» <HR«RP-C»C) / <2“bB) “*2) “TNDIF-.5«RR*S/BB* ( ( ,5“C*RR 
ll/(?“PRI““?l»RATLN 

RSTAN = RSTAN/RR 

RS2TAN = <S»»3/3) “TNSUM* ( <RR« «3-3“C“«2“RR > / (4R»BB“*3 ) ) *RATLN*C*RR* 
IS/ (3“RR*»2) - ( (ft»C“PP*»2-2“C*»3> /<6» (2«BB) *“3) I “TNOIF 
P52TAN = PS2TAN/PP 

RS3TAN = (S«“4/4) “TNOIF- ( (C“RR»*3-C*« 3*RR ) / ( 32“BB*“4 ) 1 “RATLN-TNDI F 
l«(PR*»4-e»C«“2“RR“»2*C*»4) / (64«BB#»4) -RR«S**3/ ( 12*BB1 -S“ (9“C**2“PR 
2**?-3“RR““4 1 / (48»RP*RH**3) 

PS3TAN = RS3TAN/PP 

RS4T AH = (S“«5/ (5«RR) ) « TMSUM ♦ TNO IF • 15«C«RR“*6- 1 0*C*“3“RP*“4 +C*“5*R 
1R**2 1 / ( lft0*BH*“5“RP»*3) -RATLN» ( RR««6- 1 0“C*“2«RR«“4*5“C*»4“RR““2 ) / ( 
?320“PP»»2*t'P“»5) ♦C»S*»3/ ( 15*HH»*2) *S» IC»*3-C*RR*“2) / nO“BB»*4) 

GO TO 13 

FOP CASE OF RR=TT=0,0, If 1=J 

7 TOP = C*2.“HH“S 
EOT = C-2.»PH“S 

IF (C.FQ.O.O.ANO.BP.EQ.O.O) 60 TO 13 
SBAO = -C/(?.“BH) 

SBADAB = AbS(FPAO) 

IF (SPAOAH.LT.S) GO TO 8 
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GO TO 4 F 127 

fi CONTINUE F 128 

WRITE (0.25) F 129 

130 SUL = SBAD-TOl F 130 

SLL = sea[)*T0l F 131 

CL OGB 1 = ALOG(TOP/<C.?.«HB*SLU ) F 13? 

CLOGR? = ALOG( (C.?.»HH»SUU/eOT) F 133 

CLOGR = CL0GR1 + CLOGR? F 134 

135 RTAN = - ( .5/PR) oCLOGR F 135 

RSTAN = ( .?5»C/RH»<»?) oCLOGR- ( ,5/BB) • ( 2 . »S-SLL«SUL ) F 136 

RS2TAU = -(C< M1 ?/(8.<*BH»«3) ) »CL0GR ♦(C/(4.*BB*»?))* (2 , •S“SLL*Sl)L) - ( . F 137 

12S/HH)« (SUL«SUL-SLL«SLL) F 138 

RS3TAN = (C«*3/(I6. l, BB**4))»CL0GR-{2,»S e# 3“SLL**3*SUL**3)/(6.»8B)“ F 139 

140 1 (C*C/ (R.»PH»»3) ) »(2.»S-SLL.SUL) ♦ <C/(B.«Bh«BB) t*(SUL*SUL-SLL»SLL| F 140 

PS4TAN = -(C««4/(32.»HH»»5) ) «CLUGR- ( SUL»» 4-SLL«»4 ) / ( 8 ,«BB ( ♦ (C/ ( 1 2 . F 141 
1*PR«HP) )• (2.«S«»3-SLI »»3* SUL * »3 ) -C*»2* ( SUL« »?-SLL»SLL> / 1 1 6. »BB**3 ) F 14? 
2*C««t3»(2.»S-SLL*SULl/(16.»B&<*»4) F 143 

GO TO 13 F 144 

145 9 CONTINUE F 145 

IF (ARS (TOP) .IT . 1E-9.0R.ARS (BOT) .LT.1E-9) GO TO 12 F 146 

IF (TOP. LE. 0.0) GO TO 11 F 147 

IF (POT. LE. 0.0) GO TO 11 F 148 

10 CLOGR = AI.OG ( TOP/HOT ) F 149 

150 RTAN = -(l/RB)»CLOGR F 150 

RTAN = RTAN/? F 151 

PSTAN = (.?5»C/HR«<*?)»CLOGR-S/HB F 15? 

PS2T AN = -(C»«?/(R'*PR«»3) )»CLOGR*C»S/(2.«BB»«2) F 153 

PS3T An = (C»3/(16.«HH»«4) )*CL0GR-S» # 3/(3.*BB)-S*C»»?/(4.«BB«»3) F 154 

155 RS4TAM = - (C««4/ (3?*BB<>»5> ) »CL0GR*2.«C«S*»3/ (1?»BB««2) ♦2.*S*C»»3/( F 155 

1 16i>80»<*4 I F 156 

GO TO 13 F 157 

11 TPOEiBT = TOP/POT F 158 

IF ITPPBHT ,GT ,0,0) GO TO 10 F 159 

160 IF (TOP. LT. 0.0) GO TO 14 F 160 

IF (HOT. LT. 0.0) GO TO 14 F 161 

1? CONTINUE F 162 

SAWAY = S-TOl F 163 

TOP = C«2.«e8»SAWAY F 164 

165 BOT = C-? . # BB W S AWA Y F 165 

CLOGR = AL OG ( T OP/BOT ) F 166 

RTAN = -(l/HH)*Cl.OGR F 167 

RTAN = RTAN/? F 168 


On 

On 
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RSTAU = t .?5»C/HR»«2) "CLOGR-S/HB 

F 

169 

170 


RS2TAN = -(C*«?/(H*M0«*3) )»CL0GR4C*S/I2.*88«*2) 

F 

170 



RS3TAM = (C«"3/ ( 16.*PH«»4 ) j «CL0GR-S»« 3/ ( 3 . »UB 1 -S»C* *?/ ( 4 . «BB**3 ) 

F 

171 



PS4TAN = - (C*»4/ (32*BH»«S) ) «CLOGR*2 . »C<*S« # 3/ ( l?*BRo*2) «2.*S*C»»3/ ( 

F 

172 



116*HB#«4) 

F 

173 



GO TO 13 

F 

174 

175 

13 

CONTINUE 

F 

175 


C 


F 

176 


C 

WRITE STATEUFNTS GO MERE IF NEEDED 

F 

177 


C 


F 

178 



GO TO 15 

F 

179 

180 

14 

WRITE (6,26) 

F 

180 


15 

CONTINUE 

F 

181 



GO TO 24 

F 

182 


C. 


F 

183 


C 

FOR CASE OF PR.TT NOT ZERO. BUT WITH MIDRANGE SINGULARITY 

F 

184 

105 

C 


F 

185 


16 

CONTINUF 

F 

186 



WRITE (6.27) 

F 

187 



SUE = S7EPO-TOL 

F 

188 



SLI. = S7ERO* TOL 

F 

189 

190 


SMIOl = S-0.5«A6S(S-5LL) 

F 

190 



SMID2 = -S«0.5»ABS (-S-SUL) 

F 

191 



AN6 1 = ATAN2 ( (C,2.»BR*S) , ABS (RR ♦ TT «S)) 

F 

192 



ANG2 = ATAN2 ( ( C *2 . •PB»SM I D 1 ) . ABS (RR*TT*SMID1 1 > 

F 

193 



ANG3P = AT AM2 I (C«2 • *BB*SLL) » ABS (HP*TT*SLL) 1 

F 

194 

195 


ANG3 = ATAN2 ( (C-2.*HH»S) .ARS (RR-TT»S1) 

F 

195 



ANG4 = ATAN21 (C*2.»HB»SM1D2) , ABS (RR*TT«SMID2) 1 

F 

196 



ANG5 = ATAN2( (C*2. e RB**SUL) ,AD5(RR*TT* SUL ) ) 

F 

197 



DO 23 1=1,2 

F 

198 



IF (I.FQ.2) GO TO 17 

F 

199 

200 


AA1 (1,1) = SUI «SUL 

F 

200 



AAJ (1 ,2) = SUL 

F 

201 



A A 1 ( 1 , 3 ) = A®1 ( 2 , 3 ) = A A 1 (3, 31 = 1. 

F 

202 



AA1 (2.1 ) = SM1D2»SMI02 

F 

203 



AA) (2,2) = SHIR? 

F 

204 

205 


AA1 (3,1) = S»S 

F 

205 



AA 1 (3,2) = -S 

F 

206 



A A ( 1 ) = A NO 5 

F 

207 



AA ( 2) = ANG4 

F 

20B 



AA (3) = AUG3 

F 

209 

210 


CLOGR = ALOG (AOS( (RR*TT*SUL)/(RR-TT»S))) 

F 

210 
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SUSE = SM 1 D2 F 2U 

PEIS = SUl *S F 212 

PELS? = SUL»®?-S<">? F 213 

DELS3 = SIIL»»3*S**3 F 21* 

215 DELS* = SUl *»*-s*>«* F 215 

DELSS = SUL»«s*s®«5 F 216 

0ELS6 = SUL««t-S»»6 F 217 

GO TO 18 F 218 

17 CONTINUE F 219 

220 A A 1 (1,1) = S«S F 220 

AA ltl«?) = S F 221 

A A 1 (1,3) = AA1 (2. 3)=AA1 (3,31 = 1. F 222 

AA l (2*1) = SMini«SM!01 F 223 

AAl (2,2) = SMIPI F 22* 

2?5 AAl (3,1) = SLL»SLL F 225 

AAl (3,2) = Sll F 226 

A A ( 1 ) = AND 1 F 227 

AA ( 2 ) = ANC'2 F 228 

A A ( 3 ) = ANG3P F 229 

230 FLOOR = «L06(ABS( (RR,TT*S) / (RR*TT»SLl 11 ) F 230 

SUSE = SMIOl F 231 

DELS = S-SLL F 232 

DFLS2 = S*S-SLL»SLL F 233 

DE L S3 = S»“3-SLL««3 F 23* 

235 PELS* = S«»*-SLL<>«* F 235 

PFLS5 = S«»5-SLL«»S F 236 

PEL S6 = S*>«6-SLL»*6 F 237 

IB CONTINUF F 238 

CALL SIPKO (AA1,3.AA,1,PET,IP1V0T,3,ISCALE) F 239 

2*0 Cl = A A ( 1 ) F 2*0 

C2 = A A ( 2 ) F 2*1 

C = AA (3 ) F 2*2 

COhO = IC«1T»1T-C?«RH»TT,C1«RR»RH)/(TT««3) F 2*3 

CON] = (C?«TT-C!»WR)/(TT»TT> F 2** 

2*5 CON 1 1 = 0.5«C1/TT F 2*5 

CON2 = (C«TT»TT-C2»RR»TT,C1«RR»RR) / (TT»*3) F 2*6 

CON 3 = (C»PBM r»TT-C2<RP«RR*TT*Cl»RR«*3)/ (TT*»*> F 2*7 

CON21 = C2/(2»TT)-C1«PP/(2»TT»TT) F 2*8 

CON* = (CffRq»pp»TT«TT-C2<'RR»«3»TT*Cl»PR«**|/(TT»»5) F 2*9 

250 CONS = *,ftC2»TT*n l 3-*,<>Cl«RP*TT««2 F 250 

CONS = C0N5/ (1?*TT»®*) F 251 

COMO = (-C»PP»TT»TT*C?»PR»RR»TT-C1»RR*»3)/(TT»**) F 252 


Ln 

'•J 
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C0N31 = (C«TT«»3-C?»RR»TT*»2»Cl*RR*»2*TT)/t2*TT**4) F 253 

CON32 = C 1 / ( 4 °T T I F 254 

255 CONB = 20. » ( C« T T *»4 -C2«RR » TT««3*C1 «RR« « 2«T T«« 2 ) / ( 60»TT»«5) F 255 

CONS = (C«PR<t»?®TT« , >2-C2®RR<>»3»TT*Cl*PR 0 «4)/{TT ,t »51 F 256 

FOMIO = C»PR««3»T1»*>?-C2*RR«< , 4«TT*C1»HR»*5 F 257 

CON 10 = CONI 0/ (TT“«6) F 258 

COM4 1 = <C2»TT-C1»HP1/(4*TT«»2) F 259 

260 CON42 = (-C»RR»TT»»2»C2»RR<t«2«TT-Cl«RR»»3|/(2»TT»»4) F 260 

CONA = (CflRP«K»4»TT«®?-C?»PR«»5«TT«Cl»RR»»6) / <TT**7) F 261 

CONR = (C2»TT-C1«RR)/(5»TT»«2) F 262 

CONC = <-C«RR«TT«*?*C2«RR»*2«TT-Cl«RR«*3) /(3«TT**4) F 263 

COMD = CONC*PR»»2/(TT«*2/3) F 264 

265 CONE = <C®TT««?-C?»PR*TT*C1*RR««2)/{4«TT«*3> F 265 

CONF = C]/(6«TT) F 266 

CONG = CONF •RRn*?*?/ ( TT»«2) F 267 

IF (I.EO.l) GO TO 19 F 268 

IF II.F0.2.AN0. (RR*TT*SUSE) .GT.0.1 GO TO 19 F 269 

270 PTAN = RTAN-CONO»CLOGR»COM«OELS*CONll*[>ELS2 F 270 

RSTAN = R ST AN-CON? 1 «OE.LS2 ♦CON2«OttS-CON3»CLOGR«Cl*D£LS3/ I 3*TT 1 F 271 

RS2TAN = RS?TAN-CON4®CLOGR»CON5»OELS3*CON31«OEl S2*C0N32«0ELS4*C0N6 F 272 
1 *0ELS F 273 

RS3TAN = RS3TAN-CON4 1 »DELS4 ♦C0N8*DF.L53*r0N4?*DEL5?+C0N9 # DELS-C0Nl 0 F 274 
275 1°CC. 0GR*C1»I)EI.S5/(5®TT) F 275 

RS4TAN = PS4TAN-CONA«CLOGP*CON8«OELS5*CONE»DEL54*CONF*DELS6*CONC*D F 276 
lFLS3*C,ONG*nELS?*CONn®nELS F 277 

GO TO 22 F 278 

19 CONTINUE F 279 

280 RT AN = RTAN*CONO»CLOGR*CON1«DELS*CON11»DFLS2 F 280 

RSTAN = RSTANtCON?l«nELS2*CON2«OELS-CON3»CLOGR»Cl*DFLS3/(3»TT> F 281 

RS2TAN = PS2TAN*C0N4»CL0GP*C0N5»DELS3*C0N31*I)ELS2*C0N32*0ELS4*C0N6 F 282 
1«0ELS F 283 

RS3TAN = RS3TAN*CON41»OFLS4*CON8«OFl.S3*CON42«OELS2*CON9*DELS-CON10 F 284 
285 1»C10GR*C1*0M S5/(S«TT) F 285 

PS4TAN = RS4TAN*C0NA«CL0GP*C0N8»UtLS5*C0NE«DFl.S4«C0NF*DELS6*C0NC«D F 286 
lFLS3+C0NG«0fcLS2*C0N0«0El.S F 287 

IF (I.FO.l.AMn. (RP«TT*SUSE) .LT.O.l GO TO 20 F 288 

GO TO 21 F 289 

290 20 RTAN = -RT AN F 290 

RSTAN = -RSTAN F 291 

ps?TAN = -RS2TAN F 292 

RS3T AN = -PS3TAN F 293 

RS4 T AM = -PS4TAN F 294 
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295 

21 

CONTINUE 

F 

295 


22 

CONTIMUF 

F 

296 


23 

CONTINUE 

F 

297 



GO TO 13 

F 

29B 


2« 

CONTINUE 

F 

299 

300 


RR = RPR 

F 

300 



S = SSS 

F 

301 



C = CCC 

F 

302 



RETURN 

F 

303 




F 

304 

305 

25 

FORMAT(3(IX.IlH?On ENTERED) 

F 

305 


26 

FORMAT (30X»43P0NE OF THE ENDPOINTS HAS A NEGATIVE LOG ARG) 

F 

306 


27 

FOPMAT(30x,10HP0 EMTEREO) 

F 

307 



ENO 

F 

308 


Ol 

VO 
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10 


15 


20 


25 


30 


35 
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SUBROUTINE LOGS (Stf ,U,RELN,RE5LN,PES2LN,RE53LN) 

C 

C SUPPORT IM LOGS 

C 

C CALCULATES INTEGRALS OF FORP S««N*AL0G (S*S*E»5*D) WITH 

C RESPECT TO S OVER LIMITS OF -S TO S FOR N=0.1.2«3. 

C 

C S= PANEL SEGMENT HALFWIOTH 

C F ,0 ARE CALCULATE!' IN SUBROUTINE CONCAL 

C INTEGRAL RF SI II TS ARE RELN , RE SLN . RES2LN , PES3LN 

C 

C EVALUATION OF INTFGPALS PERFORMED USING MAC5YMA ALGEBRAIC 

C MAN IPI!I ATI ON PROGRAM OF MIT PROJECT MAC 

C IF 1=J INTEGPAL EVAl.UALTEO AT APPROXIMATE ENDPOINTS » +-S AW AY 

C 

PEAL LATP,LAne,Ll.L2.L3,LA 
COMMON /TELL/ TOL.TOL2 
PFLN = 0. 

RE SLN = 0. 

RES2LN = 0. 

RFS3LM = 0. 

SS = S 

A = S*S*E*S*0 
P = S»S-E»S*0 
A = ABS(A) 

P = A0S (6 ) 

AA = ARS(A) 

ee = aps(6i 

IF (AA.LE.o.nnnnooooi > go to 6 
IF (BR.LE .0.000000001 ) GO TO 6 
1 DISC = f«F-4*D 

D I SO = SORT (APS (RISC) ) 

OIS = E»F-2.«!) 

0153 = E«*3-3.«0»t 

0154 = (E»f-4.»0) • (F»E-D) 

0IS44 = E»»4-4.»0*E»E + 2.»n'>R 

0155 = t «»6-6.«n»F <t»3«fl.»F»D«*2 
LATH = AL06(A«e> 

l ADH = ALOG ( A/B 1 

IF (AA.LE .0.000000001 ) S = SAWAY 
IF (BQ.Lt .0.000000001) S = SAWAY 
PE = S«LATEU0.t»f«LAOR 
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ESP = F*2«S G 43 

ESP = F-2»S G 44 

IF (0 ISO) 4,3.2 G 45 

? CONTINUE G 46 

LI = F-I)ISO*?«S G 47 

L? = E*UISCJ4?«S G 48 

l 3 = F-PISQ-?*S G 49 

L4 = F*PIS()-2*S G 50 

PIFFLN = ALOG(L]«L4/L2/L3> G 51 

RFLN = Rt-4»S-0.5«PISO»PIFFLN G 52 

RE SEN = n.5«LAnB«S»»?* < 0 , 25«E »l) ISC/01 SO ) «DIFFLN G 53 

RESLN = RESLN-0.?5»[iIS«LADB4E»S G 54 

RES2LN = (S«»3/3)«l ATB-0IS4/(6»0ISQt«DIFFLN4LADB»»niS3/6)-4*S«»3/9 G 55 

l-6»DIS»S/9 G 56 

RE53LN = 0.?5«S»»4»LAnB4DISb/ (8«DIS0)«DIFFLN-LA0B»(0IS44/8( 4E»S»*3 G 57 

1/640. 5«S«OIS3 G 58 

GO TO 5 G 59 

3 CLOGPT = ALOGlESP/eSM) G 60 

RELN = S»t.ATB-4»S4F'>CL0GRT G 61 

RESLN = 0.S«S« B 2<*LAPR-0,5*n I S*CLOGRT 4E*S G 62 

RES2LN = (S*«3 B LATB4PIS3»CL0GHT4 (OIS4)* ( (1/ESP>-(|/ESM) )-2»S«DIS)/ G 63 

13-4«(S»»3)/9 G 64 

RES3LN = 0.25»S»«4ttl AOH-O I S4 4 »0 . 25*CLOGRT- . 25«D IS5* 1 1 /ESP- 1/ESM) 4E G 65 

1*S«»3/640.5»S«PIS3 G 66 

GO TO 5 G 67 

4 TNRAT = ATAN? (ESP.DISO) -ATAN2 (ESH.OISQ) G 68 

RFLN = PE-4»S- (PISC/DISO) »T6PAT 6 69 

RESLN = 0 . 5 *(S«‘> 2 - 0 . 5 »OISI»LAnB 40 . 5 *E«DISC/DISO»TNRAT 4 E»S G 70 

RFS2LN = S»«3/3«LATR4 (0IS3/6) B LAOB- (OIS4/ (3»DISQ) )«TNPAT-4*S*«3/9- G 71 

lP«S B PIS/3 G 72 

RFS3LN = (0,?5 B S««4-PIS44/8)»LAOb40.?5*PIS5/DISQ*TNRAT4E*S»*3/64S* G 73 

1P1S3/2 G 74 

5 CONTINUE G 75 

GO TO 7 G 76 

6 CONTINUE G 77 

SAHAY = S-TOL? G 78 

A = SAWAY»SAWAY4E»SAWAY40 6 79 

b = SAWAY*SAWAY-F*SAWAY40 G 80 

A = ABS(A) G 81 

P = AHS(R| G 82 

GO TO 1 G 83 

7 CONTINUE G 84 



SUBROUTINE LOOS 


PAGE 


B5 


S = SS 
RETURN 
END 


G 85 
G 66 
G 87 


SUBROUTINE DRACAL 


PAGE 


1 


1 


5 


10 


15 


20 


25 


30 


35 


AO 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE ORACAL ( 1 . J , YHH , ZHH . PPP , A INT ) H 

H 

SUBROUTINE DPACAL H 

H 

TREFFT7 RLANF DRAG ANALYSIS ASSUMES PIECEWISF LINEARLY VARYING H 

SHEO VORTICITY SHEET STRENGTH H 

H 

CALCULATE INTEGRALS A THROUGH F FOR DRAG COF.F CALCULATION H 

H 

CALLS SUBROUTINES LOGS.SNTAN ,CCAL ,CONCAL H 

H 

DIMENSION AIM (6) H 

DIMENSION YHH ( 1 ) • 7HH(1), PPP(l) H 

COMMON /SEG/ SNN(Sl) H 

PI = A.«ATAN(1.) h 

S = SNEKJl H 

CALL CCAL (I.J.YHH,ZHHiPPP*S.AA,H8,DO»FF.GGtEE t AJ,AK.RRtTT.UUtWW) H 
CALL CONCAL ( A A , HH , FF , GG , S « A , B t C , D,F ,G t C J.CK *CL • CM.CN ,CO t CP 1 1 1 H 

S = SNN(I) H 

CALL LOGS |S,CJ,F,PEI.N,RESLN#RE52LN,RES3LN) H 

call SNTAN (S,C,BB,RR,TT,RTAN,RSTAN,RS2TAN,RS3TAN,RS4TAN> h 

AAAAAA = A»RTAN*8<*RSTAN*RB»PELN/4 H 

BRPRHR = 2,*(CL«RTAN*CM»RSTAN«CN*RS2TAN>*C0»PELN*CP»RESLN H 

CCCCCC = A«PSTAN*R»PS2TAN+Rfi«RESLN/A H 

DOODDIl = ?.* (CL»RSTAH*CM»PS2TAN*CN«RS3TANl ♦CO»RESLN*CP»RES?LN H 

EEEEtF = A»RSPTAN*R«PS3TAN*8MRES2LN/4 H 

FFFFFF = 2.*(CL*RS?TAMCM#PS3TAN*CN»RS4TANI«C0*RES2LN«CP»RES3LN H 

CALL LOGS (S«CH,G,RFLN,PFSUN,RES2LN,PES3LN) H 

CALL SNTAN ( S . 0 . BP , RP , T T , PTAN . RST ANt RS2T AN .RS3T AN,RSATAN) H 

AAAAAA = A AAA AA-A»RT AN-B«RSTAN-UB«REIN/4 H 

BBRFBP = HBPPPP-P.e (CL»RTAN*CM*RSTAN«CN«RS2TAN)-CO*RELN-CP»RESLN H 

CCCCCC = CCCCCC-A«RSTAN-B»RS2TAN-BB*RESLN/4 h 

oonpon = l)Unnon- 2 .» (CL»RSTAN»CM»RS2TAN+CN»RS3TAN)-CO*RESLN-CP*RES2 H 

ILN H 

EEFFtF = EEFFEF-A«RS?TAN-B»RS3TAN-BB»RES2LN/4 H 

FFFFEF = FFFFFF-P.a (CL »RS2T AN «CM*RS3T AN«CN*RS4T AN) -CO»RES2LN-CP*RE H 
1S3LN M 

S = SNN(J) H 

CALL CONCAL (nn,EF,AJ,AK,SiA,H.C,D«F,G,CJ,CK,Cl.*CM,CN,CO,CP,2) H 

S = SNM(I) H 

CALL LOGS (S,CJ,F ,RFLN , RESLN , RE S2LN , RES3LN) H 

CALL SNTAN (S.C.EE , IIU . WW , RT AN «RST AN .HS2T AN. RS3T AN .RS4TAN) H 
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Ei 6BPBB-?.*(CL*RTAN«CM»RSTAN*CN«P52TAN) -CO«tRELN-CP«RESLN 
CCCCCC -A»HSTAN-B»RS2TAN-EF »RE5LN/4 

DODPOP-?.* <CL»RSTAN«CM»RS2TAN*CN«RS3TAN) -C0»RESLN-CP»RES2 
EEEFFF-A«RS?TAN-R»PS3TAN-EE#RES2LN/4 

FFFFFF-2.*(CL»RS2TAN«CM»RS3TAN.CN»RS*TAN)-CO*RES2LN-CP»RE 


AAAAAA ♦A«RTAN*H«RS7AN*tt»RELN/4 
RFRPBR«2.* (CL»RTAN*CR«RSTAN*CN*RS2TAN) ♦CO»RELN*CP*RESLN 
CCCCCC»A«HSTAN*R*R52TAN+EE ORESLN/4 

OnnnPD*2.»(CL»RSTAN*CH*RS2TAN*CR*RS3TAN) ♦CO»RESLN*CP»RFS2 
EEEF.EE *A«R52 T AN ♦RaPS 3TAH»EE«RES2LN/4 

FFFFFF*?.* (CL«RS2TAN*CM»RS3TAN«CN*RS4TAN) *CO*PES2LN*CP*RE 
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AINT (6) 
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END 
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SUBROUTINE SIMFO ( A ,N , H , M , OE TE RM , I P I VOT ,NMAX , I SCALE ) 

I 

1 


C 


SOLUTION OF SIMULTANEOUS LINEAR EOUATIONS 

I 
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C 



I 
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DIMENSION IPIVOT(N), A(NMAX*N), B (NMAX , M ) 

I 

5 




EQUIVALENCE (IROW.JROW), ( ICOLUM, JCOLUM) , (AMAX *T .SWAP) 

I 
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C 
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INITIALIZATION 
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ISCALE = 0 

I 

10 




R 1 = 10.00*100 

I 

11 




P2 = 1.0/HI 

I 

12 




DE TERM = 1.0 

I 

13 




DO 2 J=1,N 

I 

1* 

15 


2 

I P I VOT ( J 1 = 0 

I 

15 




DO 3fl 1 = 1, N 

I 

16 
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SEARCH f OR PIVOT ELEMENT 

I 
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AMAX = 0.0 
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DO 7 J= 1 »N 

I 

21 




IF (IPTVOT IJl-l ) 3,7,3 

I 

22 
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DO 6 K=l ,N 

I 

23 




IF (IPIVOT(K)-l) 4,6,39 

I 

24 

25 


4 

IF (A8S(AMAX)-ABS(A(J,K)1) 5,6,6 

I 

25 



5 

1 ROW = J 

I 

26 




ICOLUM = K 

I 

27 




AMAX = A(J,K) 

I 

28 




CONTINUE 

I 

29 

30 
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CONTINUE 

I 

30 




IF (AMAX) 9,8,9 

I 

31 
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DETERM = 0.0 

I 

32 




ISCALE = 0 

I 

33 




60 TO 39 

I 

34 

35 


Q 

IP1V0TI ICOLUM) = 1RIV0T ( ICOLUM) *1 

I 

35 


C 



I 

36 
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INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

I 

37 


C 
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38 




IF ( I ROW- I COLl.'M ) 10,14,10 

I 

39 
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DFTERM = -DETFPM 

I 

40 




DO 11 L=1,N 

I 

41 




SWAP = A(IROVi,l ) 

I 
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A ( I ROW »L ) = AdCOLUM.L) I *3 

11 A(ICOLUM,L) = SWAP 1 44 

45 IF (M) 14,14,12 I 45 

12 00 13 L = 1,M I 46 

SWAP = H(IROW.L) I 47 

B ( I ROW * 1. 1 = R ( I COLOR , L ) I 48 

13 BdCOHIM.L) = SWAP I 49 

50 1 A PIVOT = A ( I COL IJH * ICOLUM ) I 50 

IF (PIVOT) 15,8,15 I 51 

C I 52 

C SCALE THf DETERMINANT I 53 

C I 5* 

55 15 PIVOTI = PIVOT I 55 

IF (ARS(OETERM)-Rl) 18,16,16 I 56 

16 DE TERR = 0ETEPM/R1 I 57 

ISCALF = ISCALE ♦ 1 J 5B 

IF (ABS(OETERR)-Rll 21,17,17 I 59 

60 17 DE TERM = DETERM/R1 I 60 

ISCALF = ISCALE* 1 I 61 

GO TO ?1 I 62 

18 IF (ARS(0ETERP)-R2) 19,19,21 I 63 

19 DETERM = DETEPM*R1 I 6 * 

65 ISCALF = I SC ALE - 1 J 65 

IF (ARS(0ETERM)-R2) 20,20,21 I 66 

20 DE TERM = DE TEPM«R1 I 67 

ISCALE = ISCALE -1 I 6 B 

21 IF (ARS(PIVOTl)-Rl) 24,22,22 1 69 

70 22 PIVOTI = PlVOTI/Rl I 7 0 

ISCALF = 1SCALE+1 I 7 i 

IF (ARS(PIVOTI)-Rl) 27,23,23 I 7 2 

23 PIVOTI = PIVOTI/R1 I 73 

ISCALF = 1 SCALE ♦ 1 I 74 

75 GO TO ?7 I 75 

24 IF (APS (PIVOT I ) -R2) 25,25,27 I 76 

25 PIVOTI = P1V0TI»R1 I 77 

ISCALE = I SCALE - 1 . - I 78 

IF ( ARS (PIVOTI ) -R2) 26,26,27 1 79 

RO 26 PIVOTI = PIVOTI *R1 1 80 

ISCALF = ISCALF-1 I 8 1 

27 OFTERM = DE TFRM«P I VOT I I 82 

C 183 

C D I V I OF PIVOT ROW RY PIVOT ELEMENT I 84 
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DO 29 1=1, N 

I B6 




IF (IPIVOTCU-ll 28,29,39 

I 87 



28 

A ( I COLUM «L I = A ( 1 COLUM »L I /PIVOT 

I 88 



29 

CONTINUE 

I 89 
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IF |M> 32,32,30 

1 90 
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00 31 L=1,M 

1 91 
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B(ICOLUM,L) = P(ICOLUM, D/PIVOT 

! 92 


C 
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REDUCE NON-PIVOT ROWS 

I 94 
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00 38 L 1 = 1 , M 

I 96 




IF (L 1-ICOLUP) 33,38,33 

I 97 
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T = A (L 1 , 1 COLUP ) 

I 98 




DO 35 L=1 , N 

I 99 
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IF (IPIVOTIL)-ll 34,35.39 

1 100 
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A (L 1 «L ) = A(L1.L>-A(IC0LUM,L)«T 

I 101 



35 

CONTINUE 
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IF (H| 38,38,36 

1 103 
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DO 37 1=1, M 
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37 

D ( L 1 , 1. 1 = B (LI , LI -B ( I COLUM *L) # T 
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CONTINUE 

I 106 
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RETURN 
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nonplanar wing with vertical 
fences of variable size. 
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